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Abstract 

A  three-dimensional,  unstructured,  multi-species  reacting  flow  solver  is 
developed  to  model  solid  oxide  fuel  cells  (SOFCs)  as  well  as  catalytic  reactors.  The  finite 
volume  based  solver  utilizes  density-based  method  to  solve  the  coupled  system  of 
governing  equations.  Numerical  results  for  both  SOFC  and  reactor  obtained  using  the  in- 
house  code  are  compared  with  the  experimental  results  from  the  literature  for  validation 
purposes.  Effects  of  mass  flow  rates  of  incoming  species  on  performance  of  SOFC  as 
well  as  catalytic  reactor  are  investigated.  Two  different  methods  namely  direct 
differentiation  and  discrete  adjoint  method  are  implemented  in  the  code  to  compute 
sensitivity  derivatives  for  computational  design.  The  catalytic  partial  oxidation  of 
methane  over  both  platinum  and  rhodium  catalysts  is  studied  using  the  in-house  solver. 
Eight  gas-phase  species  (CH4,  CO2,  H2O,  N2,  O2,  CO,  OH  and  H2)  are  considered  for  the 
simulation.  The  surface  chemistry  is  modeled  using  detailed  reaction  mechanisms 
including  24  heterogeneous  reactions  with  1 1  surface-adsorbed  species  for  Pt  catalyst  and 
38  heterogeneous  reactions  with  20  surface-adsorbed  species  for  Rh  catalyst.  The 
numerical  results  are  compared  with  the  experimental  data  and  good  agreement  is 
observed.  The  effects  of  the  design  variables  of  inlet  velocity,  methane/oxygen  ratio, 
catalytic  wall  temperature  and  catalyst  loading  on  the  cost  functions  representing 
methane  conversion  and  hydrogen  production  are  numerically  investigated.  The  design 
cycle  is  performed  using  two  gradeint-based  optimization  algorithms  to  improve  the 
value  of  the  implemented  cost  function  and  optimize  the  reactor  performance.  A 
capability  to  perform  thermo-mechanical  analysis  is  developed  by  coupling  the 
multispecies  solver  with  the  structures  code.  Results  obtained  using  this  capability  are 
shown  for  the  planar  SOFC.  Also,  sensitivity  derivatives  of  stress  are  computed  using  the 
direct  differentiation  method  in  the  reacting  flow  and  structures  codes.  A  novel  method  to 
perform  shape  optimization  using  CAD  based  parameters  has  been  collaboratively 
developed.  This  method  is  also  effective  in  generating  curved  elements  for  high-order 
finite  element  meshes.  A  robust  mesh  movement  algorithm  is  implemented  to  ensure  grid 
quality  during  shape  design  and  high-order  element  creation.  Examples  using  this 
technique  are  presented  for  a  tubular  SOFC,  three-dimensional  body  of  revolution  and 
NACA  4412. 

1.  Introduction 

Energy  efficiency  has  been  one  of  the  most  important  research  areas  of  recent 
times.  Improving  efficiency  of  existing  systems  as  well  as  developing  new  technologies 
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for  energy  production  and  storage  has  captured  a  lot  of  attention  lately  due  in  part  to 
reducing  reliance  on  fossil  fuels  as  well  as  their  environmental  impacts.  Systems  that  are 
capable  of  operating  with  variety  of  fuels  while  producing  energy  at  higher  efficiency  are 
very  desirable.  Solid  oxide  fuel  cells  (SOFCs)  can  be  classified  as  one  of  such 
technologies  due  to  its  capability  of  producing  efficient  energy  using  various 
hydrocarbon  fuel  types.  SOFCs  are  high  temperature  fuel  cells  capable  of  producing 
electrical  energy  at  high  efficiency.  As  the  name  suggests,  the  SOFC  contains  a  solid 
electrolyte  that  is  capable  of  transporting  oxygen  ions  from  the  cathode-electrolyte 
interface  to  the  anode-electrolyte  interface. 

Fuel  reformer  is  one  of  the  most  important  components  of  the  SOFC  system.  The 
purpose  of  the  fuel  reformer  is  to  convert  the  chemical  composition  of  primary  fuel  into 
the  species  that  systems  like  SOFC  can  be  operated  with.  Fuel  reforming  can  be  broadly 
classified  into  three  categories  including,  steam  reforming  (SR),  partial  oxidation  (POX) 
and  autothermal  reforming  (ATR).  The  reactors  used  in  reforming  process  can  have  many 
different  structures  such  as  pack  bed  and  monolith  depending  on  application  and  other 
parameters.  These  reactors  are  categorized  as  the  catalytic  reactors.  The  catalytic  reactor 
can  be  distinguished  from  the  conventional  reactor  by  considering  fundamental 
differences  between  homogeneous  (conventional)  combustion  and  catalytic  combustion. 
The  main  differences  can  be  summarized  as  [1]  following. 

•  Conventional  combustion  occurs  in  the  presence  of  a  flame,  while  catalytic 
combustion  is  a  flameless  process. 

•  Catalytic  combustion  generally  proceeds  at  a  lower  temperature  than  conventional 
combustion. 

•  Catalytic  combustion  results  in  lower  emission  of  oxides  of  nitrogen. 

•  Conventional  combustion  can  only  exist  within  well-defined  fuel  to  air  ratios. 
Catalytic  combustion  is  not  constrained  by  such  condition. 

•  Cataly  tic  combustion  can  offer  fewer  constraints  on  reactor  design. 

Numerical  simulation  techniques  have  become  matured  enough  to  be  utilized  in 
performing  design  and  optimization  in  vast  variety  of  fields.  Some  of  the  advantages  of 
the  numerical  simulations  over  the  experiments  are  the  cost  effectiveness  and  the  fact  that 
the  simulations  provide  a  wealth  of  data  that  is  difficult  or  impossible  to  obtain 
experimentally  and  can  be  used  to  perform  in-depth  analysis  of  the  system.  Numerical 
design  is  an  iterative  process  starting  with  baseline  solution,  followed  by  sensitivity 
computation  and  parameter  optimization.  Each  of  the  steps  mentioned  can  be  performed 
with  different  methods  from  high-fidelity  to  low-fidelity  as  well  as  in  one-dimension  to 
multi-dimensions.  As  expected,  high-fidelity  simulations  provide  detailed  information 
about  the  flowfield,  but  comes  with  a  burden  of  higher  computational  cost  compared  to 
the  simplified  models.  Even  though  SOFCs  are  still  in  the  developmental  stage, 
numerical  techniques  can  be  effectively  utilized  to  find  solutions  to  many  of  the  design 
hurdles  affecting  the  commercial  application  of  the  technology.  Numerical  simulations 
can  contribute  greatly  toward  the  development  of  better  designs  that  can  produce  more 
power,  increased  efficiency,  and  extended  life  expectancy  of  various  SOFCs  and 
reformers. 
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Behavior  of  SOFCs  has  been  studied  by  different  researchers  using  both 
simulation  models  as  well  as  experiments  [2-9].  Primary  focus  of  numerical  simulations 
is  on  analysis  of  fuel  cells  without  putting  much  emphasis  on  formal  optimization 
procedures  or  sensitivity  computation.  Most  of  the  design  studies  are  performed  by 
simply  changing  the  parameter  of  interest,  re-running  the  simulation,  and  comparing  the 
results  with  those  from  the  original  simulation  [2, 6, 7,9].  While  this  approach  can  be  used 
to  determine  the  effects  of  parameter  variations  on  fuel  cell  performance,  a  more  rigorous 
approach  toward  optimization  would  likely  lead  to  better  designs,  and  can  also  provide 
improved  insight  into  the  parameters  affecting  the  performance  of  the  fuel  cell.  For 
SOFC  problems,  example  cost  functions  that  can  be  used  for  improving  performance 
include  minimizing  temperature  variations,  obtaining  equal  distribution  of  fuel  in  each  of 
the  channels,  minimizing  stress  inside  different  components  or  maximizing  power. 

Design  variables  may  be  related  to  the  shape/size  of  the  fuel  channels,  electrodes, 
electrolyte,  and  interconnect,  but  may  also  be  coupled  to  the  stoichiometric  composition 
of  fuel  or  material  properties  such  as  the  porosity  or  tortuosity  of  the  electrodes.  In 
references  [10]  and  [11],  optimization  algorithms  have  been  used  to  improve  the 
performance  of  a  polymer-electrolyte-membrane  fuel  cell  (PEM)  using  four  design 
variables,  where  the  sensitivity  derivatives  used  for  the  optimization  algorithm  have  been 
obtained  using  a  finite-difference  approach.  While  finite  differences  are  often  a  viable 
means  for  computing  sensitivity  derivatives,  this  method  can  be  computationally 
restrictive  when  a  sufficiently  large  number  of  design  variables  are  present.  In  addition, 
accurate  derivatives  can  sometimes  be  difficult  to  obtain  using  finite  differences  because 
of  subtractive  cancellation  errors  [12],  which  occur  when  the  function  evaluations  in  the 
numerator  become  computationally  indistinguishable  [13]  when  very  small  perturbations 
are  used.  By  using  a  direct  differentiation  or  discrete  adjoint  method  [14-25],  sensitivity 
derivatives  that  are  consistent  with  the  flow  solver  may  be  obtained  for  use  in  a  design 
optimization  environment.  The  code  utilized  in  this  report  is  capable  of  computing 
sensitivity  derivatives  of  a  cost  function  with  respect  to  desired  parameters  using  both 
direct  differentiation  and  discrete  adjoint  methods.  In  different  design  studies  performed 
at  the  SimCenter,  sensitivity  derivatives  were  also  utilized  in  a  design  environment  to 
optimize  cost  functions  representing  uniform  fuel  distribution  [26,27],  temperature 
distribution  [28]  and  cell  voltage  [27]. 

Due  to  the  presence  of  various  thermal  mechanisms  inside  SOFC,  temperature 
distribution  exhibits  non-uniformity  throughout  the  domain.  This  coupled  with  the 
mismatch  in  coefficients  of  thermal  expansion  of  different  SOFC  components  makes 
stress  analysis  an  interesting  avenue  to  explore.  Recently,  few  studies  have  been 
performed  [29-35]  to  analyze  stress  components  inside  different  components  of  SOFC. 
Lin  et  al.  [29]  analyzed  effects  of  clamping  load  on  the  thermal  stress  distribution  in  a 
planar  SOFC.  Five  different  compressible  loads  were  applied  to  investigate  effects  on 
stress  distribution.  Gulfam  et  al.  [32]  analyzed  thermal  stress  inside  PEN  region  of  SOFC 
for  co-flow,  counter- flow  and  cross-flow  configurations  using  commercial  software, 
ABAQUS  [36].  Maximum  stress  inside  anode  layer  was  found  at  high- temperature 
regions  located  on  the  anode-electrolyte  interface  for  all  configurations.  Jiang  et  al.  [33] 
performed  thermal  stress  analysis  of  SOFC  with  the  bonded  compliant  seal  design.  Stress 
analysis  performed  using  commercial  software,  FLUENT  [37]  and  ANSYS  [38] 
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investigated  effects  of  temperature  non-uniformity  and  cell  voltage  on  thermal  stress 
distribution.  Weil  and  Koeppel  [34]  also  analyzed  effects  of  different  seal  designs  on 
stress-strain  state.  They  used  glass-ceramic,  brazed  joints  and  foil-based  seals  in  this 
study.  Chiang  et  al.  [35]  analyzed  effects  of  anode  porosity  on  thermal  stress  in  anode- 
supported  SOFC  using  commercial  software,  STAR-CD  [39]  and  MARC  [40],  The  study 
indicated  presence  of  higher  principal  stress  at  low  cell  voltages  due  to  high  local  current 
density  and  steep  temperature  gradients.  None  of  these  studies  attempted  to  compute 
sensitivity  derivatives  of  a  cost  function  involving  thermal  stress  components  using 
formal  procedures  such  as  direct  differentiation  or  discrete  adjoint  methods.  Such 
capability  is  required  if  optimization  of  principal  stress  components  or  strain  rates  with 
respect  to  geometrical  or  material  parameters  is  desired.  One  of  the  cases  described  in 
this  report  demonstrates  such  capability,  where  stress  sensitivities  are  computed  with 
respect  to  the  cathode  porosity  using  direct  differentiation  method.  Implementation  of  this 
method  is  not  trivial,  as  it  requires  computation  of  sensitivity  derivatives  of  flowfield 
variables  in  the  multispecies  Navier-Stokes  code  and  coupling  them  with  the  structures 
solver  to  compute  stress  sensitivities  with  respect  to  the  design  parameters. 

The  monolith  or  honeycomb  reactor  is  a  commonly  used  configuration  in  catalytic 
combustion.  It  consists  of  a  number  of  parallel  passageways  through  which  the  gas  flows, 
with  the  catalyst  being  coated  on  its  walls.  Catalytic  monolithic  reactors  are  generally 
characterized  by  the  complex  interaction  of  various  physical  and  chemical  processes  that 
include  transport  of  momentum,  energy,  and  chemical  species.  Modeling  of  monolithic 
reactors  can  be  broadly  divided  into  two  categories:  single-channel  modeling  that 
considers  one  channel  of  the  monolith  and  full-scale  modeling  that  considers  the  whole 
reactor  comprised  of  several  hundred  channels  [41].  Single-channel  models  can  be  one¬ 
dimensional,  two-dimensional  or  three-dimensional.  One-dimensional  (ID)  models 
ignore  radial  and  angular  gradients  in  temperature,  concentration,  and  velocity,  and 
consider  only  axial  variations.  These  models  with  lumped  heat  and  mass  coefficients 
were  widely  used  because  of  their  simplicity  and  easy  implementation.  A  simplifying 
assumption  of  plug  flow  is  frequently  made  in  ID  models.  Inside  the  monolith  channel, 
the  catalytic  reactions  occur  in  the  washcoat  on  the  channel  wall.  There  are  two  choices 
for  incorporating  these  catalyst  reactions  in  numerical  model. 

•  The  pseudo-homogeneous  model:  The  wall  temperature  and  concentrations  are 
assumed  to  be  the  same  as  fluid,  and  the  reaction  rate  is  incorporated  directly  into 
the  conservation  equations. 

•  The  heterogeneous  model:  The  gas-solid  interface  at  the  wall  is  assumed  to  be 
discontinuous  and  separate  mole  and  energy  balance  equations  are  solved  for  the 
solid.  These  equations  are  coupled  to  the  fluid  equations  through  mass  and  heat 
transfer  coefficients.  Catalytic  reactor  results  presented  in  this  report  utilizes 
heterogeneous  model  for  surface  chemistry. 

Catalytic  combustion  of  hydrogen  was  investigated  by  Cerkanowicz  et  al.  [42] 
with  simplified  chemistry  and  by  Kramer  et  al.  [43]  with  detailed  kinetics.  Multi¬ 
dimensional  models  of  catalytic  reactors  can  be  developed  based  on  either  boundary- 
layer  equations  or  full  Navier-Stokes  equations.  In  boundary -layer  approximation,  axial 
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(streamwise)  diffusive  transport  is  neglected,  but  detailed  transport  to  and  from  the 
channel  walls  is  retained.  Raja  et  al.  [44]  investigated  the  efficiency  and  validity  range  of 
the  Navier-Stokes,  boundary  layer  and  plug-flow  models  in  a  catalytic  monolithic 
channel.  Their  research  showed  that  the  boundary-layer  models  provide  accurate  results 
with  relatively  low  computational  cost  [44],  Deutschmann  et  al.  [45]  and  Dogwiler  et  al. 
[46]  used  two-dimensional  Navier-Stokes  models  with  detailed  heterogeneous  and 
homogeneous  chemistry  for  simulation  of  catalytic  combustion.  Catalytic  combustion  of 
methane-air  was  studied  by  Markatou  et  al.  [47]  using  two-dimensional  boundary  layer 
model.  Kumar  [41]  developed  a  coupled  implicit  solver  to  solve  species  conservation 
equations  and  investigated  the  flowfield  in  a  full-scale  3D  catalytic  converter.  Catalytic 
combustion  of  iso-octane  over  rhodium  catalysts  was  studied  by  Hartmann  et  al.  [48]. 
They  used  detailed  surface  chemistry  including  17  surface  species  and  58  surface 
reactions  for  this  simulation.  Maestri  and  Cuoci  [49]  used  OpenFOAM  to  simulate 
heterogeneous  catalytic  systems  in  three-dimensions  with  detailed  kinetics  schemes. 
Catalytic  partial  oxidation  (CPOX)  of  methane  over  a  honeycomb  reactor  was 
numerically  studied  by  Hettel  et  al.  [50].  They  coupled  OpenFOAM  [51]  and  DETCHEM 
[52]  to  model  a  large-scale  COPX  reactor.  Minh  [53]  developed  numerical  methods  for 
the  simulation  and  optimization  of  complex  processes  in  catalytic  monoliths  for  two 
practical  applications:  catalytic  combustion  of  methane  and  conversion  of  ethane  to 
ethylene.  In  this  study,  [54]  optimization  of  the  oxidative  dehydrogenation  of  ethane  to 
ethylene  over  platinum  was  investigated  using  a  two-dimensional  model  to  simulate  a 
single  monolith  channel.  In  the  work  presented  in  this  report,  an  in-house  three- 
dimensional  multispecies  solver  is  coupled  with  Cantera  [55]  to  solve  for  the 
heterogeneous  reactions  present  in  catalytic  reactors.  The  solver  is  also  coupled  with  the 
DAKOTA  [56]  toolkit  to  perform  optimization  of  the  reactor. 

Computational  fluid  dynamics  (CFD)  methods  are  generally  classified  as  two 
distinct  families  of  schemes:  pressure-based  and  density-based  methods.  The  pressure- 
based  algorithm  solves  the  momentum  and  pressure  correction  equations  separately.  The 
density-based  solver  solves  the  governing  equations  of  continuity,  momentum,  energy 
and  species  transport  simultaneously.  In  density-based  approach,  velocity  field  is 
obtained  from  the  momentum  equations  and  the  continuity  equation  is  used  to  obtain  the 
density  field.  Pressure  field  is  determined  from  the  equation  of  state  using  computed 
flowfield  variables.  In  pressure-based  methods,  since  there  is  no  independent  equation  for 
pressure,  a  special  treatment  is  required  in  order  to  achieve  velocity-pressure  coupling 
and  enforcing  mass  conservation.  Originally,  pressure-based  approach  was  developed  for 
low-speed  incompressible  flows,  while  density-based  approach  was  mainly  used  for  high¬ 
speed  compressible  flows.  However,  this  distinction  has  been  blurred  in  recent  times  as 
both  methods  have  been  extended  and  reformulated  to  solve  for  a  wide  range  of  flow 
conditions  beyond  their  original  intent.  As  majority  of  work  involving  simulation  of  the 
catalytic  combustion  uses  pressure-based  schemes,  relatively  less  research  has  been 
performed  in  this  field  using  fully  coupled  density-based  methods.  In  the  study  presented 
in  this  report,  potential  of  using  the  density-based  approach  for  solving  chemically 
reacting  flow  inside  a  catalytic  reactor  and  SOFC  is  investigated.  Since  all  governing 
equations  including  species,  momentum  and  energy  are  solved  simultaneously,  very 
accurate  solution  is  obtained.  One  of  the  drawbacks  of  the  density-based  methods  is  that 
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the  system  of  equations  becomes  very  stiff  at  low  velocity.  This  problem  is  usually  fixed 
by  using  appropriate  preconditioners. 

An  in-house  three-dimensional,  multispecies  solver  is  developed  to  obtain 
different  simulation  results  described  in  this  report.  Two  validation  cases  are  included  in 
this  report.  The  first  validation  case  compares  polarization  curve  for  a  planar  SOFC  with 
the  experimental  results  [57].  The  second  validation  case  involves  catalytic  oxidation  of 
methane.  In  the  reacting  flow,  mass  flow  rate  of  incoming  species  is  considered  an 
important  parameter.  Effect  of  different  mass  flow  rates  of  incoming  species  on 
performance  of  SOFC  as  well  as  catalytic  reactor  is  investigated.  Another  goal  is 
implementation  of  formal  methods  to  accurately  and  efficiently  compute  sensitivity 
derivatives.  Sensitivity  derivatives  for  different  cost  functions  involving  SOFC  and 
catalytic  reactor  are  computed  using  discrete  adjoint  method  and  direct  differentiation 
and  compared  w  ith  each  other  for  validation  purposes.  Effects  of  different  design 
parameters  on  reactor  performance  are  also  investigated.  The  catalytic  reactor  is 
numerically  optimized  using  two  different  gradient-based  algorithms.  One  of  the  goals  of 
this  project  is  development  of  fluid-structure  interaction  capability  to  analyze  thermal 
stress  distribution  and  stress  sensitivities.  Application  of  this  capability  is  show'n  for 
different  components  of  SOFC  in  this  report.  Finally,  results  are  shown  for  obtaining 
surface  derivatives  with  respect  to  CAD  based  design  parameters.  This  method  has  been 
developed  in  collaboration  with  the  Aerospace  Computational  Design  Laboratory  of 
Massachusetts  Institute  of  Technology  (MIT).  The  surface  derivatives  obtained  using  this 
method  are  fed  into  a  finite  element  mesh  movement  solver  to  modify  the  shape  of  the 
tubular  SOFC.  Application  of  this  capability  is  also  shown  to  accurately  generate  high- 
order  grids  when  underlying  CAD  model  has  curvatures. 

2.  Governing  Equations  and  Computational  Methodology 

2.1  Multispecies  Navier-Stokes  and  Electric  Potential  Equations 

The  three-dimensional  code  utilized  to  perform  various  simulations  in  this  project 
solves  multi-species  Navier-Stokes  equations.  For  SOFC  simulations,  this  model  is 
combined  with  an  electric  potential  equation  that  governs  the  distribution  of  electric 
potential  and  current  density  in  the  field.  To  solve  for  the  surface  chemistry'  in  reforming 
simulations,  an  interface  with  Cantera  is  developed.  For  reforming  simulations,  solution 
of  electric  potential  is  not  required,  thus  it  is  ignored.  The  three-dimensional  model 
accounts  for  all  components  of  the  SOFC,  including  the  anode,  cathode,  electrolyte, 
interconnects,  and  the  fuel  and  air  channels  as  well  as  different  components  of  a  reformer. 
Note  that  the  model  is  not  limited  to  any  particular  type  of  SOFC,  i.e.  planar  as  well  as 
tubular  type  SOFC  can  be  simulated  using  this  model.  The  code  is  general  enough  to 
solve  for  different  kind  of  reacting  flows  as  well  as  surface  chemistry  reactions. 

The  governing  equations  for  mass,  momentum  and  energy  conservation  are 
solved  simultaneously  with  the  equation  governing  the  electric  potential  in  the  numerical 
model.  The  system  of  equations  utilized  in  the  SOFC  model  is  given  by  equations  (1)  - 
(6),  which  represent  the  conservation  statements  for  the  species  concentrations, 
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momentum  (x,  y  and  z),  energy  and  current,  respectively.  Equation  (6)  is  ignored  in 
catalytic  reactor  simulations. 

Equations  (1)  -  (5)  are  modified  Navier-Stokes  equations  valid  for  both  porous 
and  fluid  regions.  Detailed  discussion  on  flux  formulation  for  these  equations  can  be 
found  in  previous  work  [26].  Equation  (6)  represents  the  electric  potential  equation.  As 
solid  regions  (interconnect  and  electrolyte)  are  considered  zero-velocity  regions,  only 
energy  and  electric  potential  equations  are  solved  inside  them.  Electric/ionic 
conductivity,  <7,  in  equation  (6)  is  a  strong  function  of  the  temperature.  Expressions 
describing  the  relationships  between  the  electric/ionic  resistivity  (reciprocal  of 
conductivity)  and  the  temperature  for  various  components  of  SOFC  are  presented  in 
Table  1  [57,58]  along  with  thermal  conductivities  and  other  material  properties  of 
different  components  of  the  SOFC. 


( dPt'  +v.(£py)+v.(j>5. 
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d(spu)  „  .  -  dP  _  .  .  £2ufl 

J+y<£puV)=  £  +V«(er )  ^ 

at  ax  B 

(2) 

3(epv)  -  dP  £2vp 

3,  +^.(epvn=  e^+V.(ery)  B 

(3) 

d(£pw)  -  dP  £2wp 

+  V.(epwF)=  £  +  V*(£T  ) 

at  az  B 

(4) 

d(<eEt) +  V'(£(E'  +P)V)+V*(£j.H.)=  V.(eprF)  V*^e//  +  Vp.(oVp) 

(5) 

V*(oV0)  =  0 

(6) 

As  presented,  equation  (6)  is  an  elliptic  equation  contrary  to  the  rest  of  the 
governing  equations,  which  are  hyperbolic-parabolic  equations.  Equation  (6)  is  solved  in 
the  entire  domain  except  for  the  fuel  and  air  channels,  which  are  pure  fluid  regions. 

Several  transport  processes  take  place  at  the  anode-electrolyte  and  the  cathode- 
electrolyte  interfaces  that  strongly  affect  the  overall  behavior  of  the  SOFC.  The 
electrochemical  reactions  taking  place  at  the  cathode-electrolyte  and  anode-electrolyte 
interfaces  can  be  described  by  equations  (7)  and  (8),  respectively. 

0.5O2  +2e~  — »  02~  (7) 

H2+02~  H20  +  2e  (8) 

The  effect  of  the  aforementioned  electrochemical  reactions  is  modeled  by 
applying  mass  flux  conditions  at  the  cathode-electrolyte  (equation  (9))  and  anode- 
electrolyte  (equations  (10)-(1 1))  interfaces  using  Faraday’s  law. 
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Jn  = — —  M 

°1  4  f  °2 

(9) 

(10) 

J»,o  =  2/7  M»,0 

(11) 

In  above  equations,  /  is  the  local  current  density  and  F  is  Faraday’s  constant.  A 
negative  sign  implies  that  the  flux  is  leaving  the  interface. 


Table  1.  Material  properties  of  various  components  of  SOFC  [57,58] 

Electric  resistivity  of  anode  (Qm) 

2.98xl0“5  exp(-1392/T) 

Electric  resistivity  of  cathode  (Q m) 

8.1  lxlO-5  exp(600/r) 

Electric  resistivity  of  interconnect  (Qm) 

6.41xl0“8 

Ionic  resistivity  of  electrolyte  (Qm) 

2.94xl0-5  exp(10350/ J) 

Thermal  conductivity  of  anode  (W  ) 

6.23 

Thermal  conductivity  of  cathode  (i W  m~xK~x ) 

9.6 

Thermal  conductivity  of  interconnect  (W  m~'K~' ) 

9.6 

Thermal  conductivity  of  electrolyte  (w  ) 

2.7 

Porosity  of  anode 

0.38 

Porosity  of  cathode 

0.5 

Tortuosity  of  anode 

1.5 

Tortuosity  of  cathode 

1.5 

Permeability  of  anode  ( m 2) 

1.0x10“'° 

Permeability  of  cathode  ( m 2) 

1.0x10“'° 

Pore  diameter  of  anode  (m) 

2.0x10^ 

Pore  diameter  of  cathode  ( m ) 

2.0x10^ 

To  account  for  the  heat  generated  due  to  electrochemistry,  heat  flux  proportional 
to  the  entropy  change  associated  with  the  electrochemical  reaction  is  applied  at  the 
anode-electrolyte  and  cathode-electrolyte  interfaces.  This  heat  flux  is  proportional  to  the 

molar  formation  rate,  (/  /  nfj ,  where  ne  is  the  number  of  electrons  participating  in  the 

electrochemical  reaction. 

In  addition  to  the  electrochemical  reactions,  two  chemical  reactions,  methane 
reforming  (12.1)  and  water  gas  shift  (12.2)  reactions  are  also  available  when  methane  is 
used  as  a  fuel  and  internal  reforming  is  allowed. 
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(12.1) 


CH  +  H  O  <: 


w. 


_ l  CO  +  3H. 

kb.  2 


CO  +  H.O  < - >  ca  +  H ,  (12.2) 

2  2  2  v 

Reaction  rates  for  various  species  have  been  computed  using  the  following 
equations. 

Rater  =  kfrpCH4pH20  -  kbrpcop3H2  (13.1) 

Rule  s  —  kfsPcoPmo  ~  kb sp  QQ2p  hi  (13.2) 


Subscripts  “  r  ”  and  “  s  ”  stand  for  reforming  and  shift  reactions,  respectively. 
Reaction  rate  constants,  kf  and  kb ,  are  computed  using  the  methodology  outlined  in 
reference  [28]. 

The  voltage  output  of  the  SOFC  strongly  depends  on  several  irreversibilities  or 
losses  encountered  in  the  flowfield  including  activation  polarization,  concentration 
polarization  and  Ohmic  polarization.  Noren  and  Hoffman  [59]  have  provided  extensive 
discussion  on  accurately  modeling  the  activation  polarization.  The  SOFC  model  used  in 
this  work  employs  the  Butler-Volmer  equation  to  compute  activation  polarization  [59], 

The  Butler-Volmer  equation  can  be  written  as. 


f  TP  \ 

n  F 

f 

i=io 

exp 

a^—ri  , 

R  j,  lac, 

V  U  J 

-exp 

V 

n  F 
(1  -a)-s—T] 
RT  1 


(14) 


The  activation  polarization  is  denoted  by  TJacl . 

a  is  the  charge  transfer  coefficient  and  assumed  to  be  0.5  in  the  current  work. 
ne  represents  the  number  of  electrons  involved  in  the  electrochemical  reaction,  which  is 
2  (equations  (7)  and  (8))  in  the  current  simulation. 

iQ  is  the  exchange  current  density  and  is  computed  using  equations  (15)  and  (16)  for  the 
anode  and  cathode  [60],  respectively. 
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’  RT  , 
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(15) 


(16) 


Various  constants  in  the  above  equations  are  given  in  Table  2  [60].  Once  the 
values  of  a  and  ne  are  inserted  in  equation  (14),  the  activation  polarization  can  be 
computed  using  the  following  expression. 
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(17) 
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Ohmic  polarization  is  a  direct  consequence  of  the  resistance  offered  to  the  flow  of 
electrons/ions  inside  various  components  of  the  SOFC.  Voltage  drop  due  to  Ohmic 
resistance  is  directly  proportional  to  the  current  and  the  resistance.  The  effect  of  Ohmic 
polarization  on  the  voltage  loss  is  directly  included  in  the  potential  equation,  equation  (6), 
through  the  electric  conductivity,  <j ,  which  is  the  reciprocal  of  the  electric  resistivity. 


Table  2.  Constants  used  to  compute  activation  polarization  [60] 

a 

0.5 

n 

t 

2 

Ca(A  rn2) 

5.5x10s 

Ze(A  m~2) 

7.0x10s 

Enr,SJ  tonoP') 

oo 

o 

i—H 

X 

o 

H 

Eari  (J  kmor') 

1.2x10s 

Prpf{N  rn'2) 

101325 

Concentration  polarization  is  caused  by  reductions  in  the  concentrations  of  the 
reacting  species  at  the  interface  between  the  electrodes  and  the  electrolyte.  The  effect  of 
the  reduction  in  concentrations  can  be  seen  from  the  well-known  Nernst  potential 
equation,  given  by  equation  (18).  Also,  exchange  current  densities  at  the  anode- 
electrolyte  interface  and  the  cathode-electrolyte  interface,  represented  by  equations  (15) 
and  (16),  respectively,  are  strongly  affected  by  the  concentration  polarization. 

Equation  (18)  computes  the  electromotive  force  (EMF)  or  electric  potential  under 
reversible  conditions,  i.e.  in  the  absence  of  activation,  Ohmic  or  any  other  losses. 


EMF  =  EMF0  +  —  In 
IF 


(  p  p°s\ 

hSo, 


V 


hj:> 


-  P 

,  where  P.  = — — 

'  P  , 

ref 


(18) 


The  electromotive  force  at  standard  pressure,  EMF° ,  is  computed  using 
polynomial  thermodynamic  relationship  between  Gibbs  free  energy  and  temperature  of 
different  species  participating  in  the  electrochemical  reactions.  The  value  of  Pr .  is  taken 

as  one  atmosphere  in  the  above  equation. 


The  electrochemical  reaction  reduces  the  concentration  of  the  reactants  and 
increases  the  concentration  of  the  products  at  the  electrode-electrolyte  interface.  Thus, 
the  partial  pressures  of  the  reactants  and  products  are  affected  in  the  same  manner.  This 
will  reduce  the  value  of  the  second  term  on  the  right-hand  side  of  the  equation  (18) 
thereby  affecting  the  EMF  of  the  cell  adversely.  Concentration  polarization  strongly 
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depends  on  the  material  properties  of  the  electrodes  that  are  responsible  for  the  transport 
(diffusion  and  convection)  of  the  reactants  and  products,  to  and  from  the  electrode¬ 
electrolyte  interface. 

2.2  Surface  Chemistry 

The  heterogeneous  and  homogeneous  chemical  reaction  mechanisms  are  key 
components  of  the  reacting  flow  modeling.  The  mechanism  of  heterogeneously 
catalyzed  gas-phase  reactions  can  be  described  by  the  sequence  of  elementary 
reaction  steps  including  adsorption,  surface  diffusion,  chemical  transformations  of 
adsorbed  species,  and  desorption.  Several  modeling  approaches  are  available  to 
compute  the  reaction  rates  of  heterogeneous  reactions.  Different  approaches  such  as, 
Ab-initio  calculation,  density  function  theory  (DFT)  and  kinetic  Monte  Carlo 
modeling  are  used  to  include  the  molecular  aspects  of  heterogeneous  catalysis.  In 
power-law  kinetic  approach,  rate  of  the  catalytic  reaction  is  calculated  by  fitting 
empirical  equations  to  the  experimental  data.  In  the  last  two  decades,  the  mean-field 
approximation  (MF)  has  been  used  as  a  work-around  in  order  to  overcome  the  much 
simpler  Langmuir-Hinshelwood  or  even  power-law  approaches  and  to  include  some 
of  the  elementary  aspects  of  catalysis  into  models  suitable  for  numerical  simulation 
of  catalytic  reactors  [61].  In  the  mean-field  approximation,  the  rate  equations  similar 
to  homogeneous  reactions  are  used  to  model  heterogeneous  reactions. 

In  the  MF  model,  the  source  term,  s(-,  of  gas-phase  species  due  to 
adsorption/desorption  and  surface  species  (adsorbed  species)  are  given  by, 


y^s 

£*k=l 


^kfku  jt:Ns 


(i  =  l,...,Ng  +  Ns) 


(19) 


Where  Ks  is  the  number  of  elementary  surface  reactions  (including 
adsorption  and  desorption),  Ns  is  the  number  of  species  adsorbed  and  Ng  is  the 
number  of  gaseous  species.  The  heterogeneous  flux  on  the  surface  is  obtained  by 
equation  (20). 


Fluxhet  =  MW^  (20) 

Since  the  catalyst  is  dispersed  as  small  particles  in  the  reactor  support,  the 
active  catalyst  area  is  usually  much  greater  than  the  geometric  surface  area.  The  ratio 
of  catalyst  and  geometric  area  is  defined  as, 

r?  _  ^catalyst  /-ii\ 

r  cat /geo  ~ 

" geometric 

To  account  for  pore  diffusion  within  the  catalyst  coating  layer,  the 
effectiveness  factor,  T],  is  defined.  Fcat/geo  and  t]  are  experimentally  determined. 
Therefore,  the  heterogeneous  flux  formula  can  be  modified  as, 


FluXfoet  —  Fcat/geor\MWiSi 


(22) 
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The  temperature  dependence  of  the  rate  coefficients  in  equation  (19)  is 
described  by  a  modified  Arrhenius  expression  shown  in  equation  (23). 

kf„  =  AkTPkexp  [^]  nfi,  e^exp  [^i]  (23) 

In  equation  (23),  0j  is  surface  coverage  for  species  “i”.  For  some  simple 
surface  reaction  mechanisms,  it  is  convenient  to  specify  the  surface  reaction  rate 
constant  in  terms  of  a  “sticking  coefficient”  (probability)  rather  than  the  actual 
reaction  rate.  This  approach  is  only  allowed  when  there  is  exactly  one  gas-phase 
species  reacting  with  the  surface. 


RT 


27 iM 


(24) 


In  equation  (24),  sf  is  the  initial  (uncovered  surface)  sticking  coefficient;  t  is 
sum  of  surface  reactants’  stoichiometric  coefficients  and  T  is  the  surface  site  density 
(mol/m:).  Using  equation  (23),  the  equation  (19)  can  be  rewritten  as. 


*t  =  2&1  vik  (AkT**exp  nfta  e^exp  [^i])  Yf£H'[Xj]9ji 

(i  =  1,  ...,Ng  +  Ns) 


(25) 


The  surface  molar  concentration  of  a  species  is  obtained  by  [62], 


[Xt]  =  ®iT  (i  =  Ng  +  1 . Ng  +  Ns) 

From  equation  above,  S;  = 


(26) 

(27) 


6 ©j  _  S( 

st  ~  r 


(28) 


The  surface  site  densities  are  of  the  order  of  10-9  mol/cm2  (approximately 
1015  adsorption  sites  per  cm2)  [63].  Equation  (28)  assumes  that  the  total  surface  site 
density  T  is  constant.  The  equation  above  is  used  for  a  transient  simulation.  In  a 
steady-state  calculation,  surface  species  concentrations  (or  site  fractions)  remain 
constant  with  time  [62],  which  gives, 


st  =  0  (i  =  Ng  +  1 . Ng  +  Ns)  (29) 

At  steady  state,  surface  species  concentrations  have  to  adjust  themselves  to  be 
consistent  with  the  adjacent  gas-phase  species  concentrations  such  that  the  condition 
Si  =  0  (equation  (29))  is  satisfied.  The  solution  of  equations  (25)  and  (29)  provide  the 
surface  coverages  and  the  surface  molar  concentrations.  Once  these  values  are  obtained, 
the  source  terms  for  the  governing  equations  can  be  computed.  The  system  of  equations 
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generated  by  equations  (25)  and  (29)  is  considered  to  be  extremely  stiff.  C-language 
version  of  VODE,  CVODE,  which  is  included  in  the  Sundials  package,  is  chosen  for 
solving  these  stiff  equations.  For  coupling  the  flow  solver  with  CVODE,  an  interface 
based  on  Cantera  [55]  is  used.  The  structure  of  this  interface  is  illustrated  in  figure  1.  A 
Cantera  input  file  is  written  based  on  the  application  including  definition  of  gas  and 
surface  phases  and  detailed  chemical  reactions.  This  input  file  is  read  and  used  to  create 
and  allocate  the  Cantera  gas,  surface  and  interface  objects  at  the  beginning  of  the 
simulation.  During  simulation,  flow  solver  gives  the  gas  phase  information  including 
temperature,  pressure  and  mole  fractions  of  the  species  to  Cantera.  Then,  Cantera  sets  the 
required  parameters  and  sends  to  the  CVODE  solver.  The  surface  coverages  and  reaction 
rates  are  computed  and  communicated  back  to  the  flow  solver  to  use  as  the  chemical 
source  terms. 


Open  and  read  input  file 
(Define  gas  and  surface  species,  reactions) 


Gas  temperature,  pressure 
and  mole  fraction  for  the 
catalytic  wall  boundary 
cell 


Flow  solver 


Create  the  gas  phase  object 

- r — ~ 

Create  the  surface  phase  obieet 


I 

>bjei 

gas 

I 


Create  the  interface  object  (interface  between 
surface-gas  phases) 


Get  the  gas  information 

T 


Set  Temperature,  pressure  and  concentration 
for  the  gas  and  surface  phases 

1 11 '  '  "  " . . 

Solve  stiff  equations  using  CVODE  and 
comnutc  surface  converges 


i 


Compute  chemical  reaction  rates 


Chemical  reaction  rates 


Figure  1.  Data  exchange  between  the  solver  and  Cantera  through  an  interface 


2.3  Boundary  and  Interface  Conditions 

To  illustrate  various  boundary’  conditions  implemented  in  the  solver,  an  example 
of  a  planar  SOFC  is  shown  in  figure  2.  The  geometry  includes  all  relevant  SOFC 
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components  including  air/fuel  channels,  anode,  cathode,  electrolyte  and  interconnects. 
No-slip,  adiabatic  wall  boundary  conditions  are  applied  at  the  interfaces  between  the 
electrodes  and  the  interconnect,  as  well  as  the  side,  top  and  bottom  walls.  Same  boundary 
conditions  are  applied  at  the  walls  of  reformers.  No-slip,  isothermal  wall  as  well  as 
catalytic  wall  boundary  conditions  are  also  implemented  and  utilized  in  catalytic  reaction 
cases.  Solid  regions  (interconnect  and  electrolyte)  are  considered  zero-velocity  regions 
and  the  only  variables  computed  inside  these  regions  are  temperature  and  electric 
potential.  A  fixed  potential  (0  =  0)  boundary  condition  is  applied  at  the  bottom  wall, 

whereas  the  top  wall  is  treated  by  specifying  average  current  density  (z  =  iappUed )  in  SOFC 
simulations. 

Inflow  boundary  conditions  with  specified  mass  flow  rate  and  species  mole 
fractions  are  applied  at  both  fuel  and  air  channel  inlets.  Specified  back  pressure  outflow 
conditions  are  applied  at  both  air  and  fuel  channel  outlets.  Same  inflow/outflow  boundary 
conditions  are  applied  in  the  reforming  case.  The  code  is  capable  of  handling  fluid, 
porous  and  solid  regions  as  well  as  interfaces  between  them.  Specific  properties  of 
different  regions  as  well  as  surfaces  are  provided  through  input  files. 


(a).  Front  view 


applied 


Current  =  0 


j,  j.l  l>  l;«  M  I  If  I  I  MM  I  H  I  1 


applied 


• 

w  . 

• 

- tn-  air 

Current  =  0 

^  Out-air 

(b).  Side  View 


Figure  2.  Boundary  and  interface  conditions  (I  -  Interconnect,  C  -  Cathode,  AC  -  Air 
Channel,  E  -  Electrolyte,  FC  -  Fuel  Channel,  A  -  Anode) 
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2.4  Structures 

As  mentioned  earlier,  capability  to  perform  thermo-mechanical  stress  analysis  has 
been  developed  to  analyze  stress  inside  different  components  of  the  SOFC  and  reformers. 
Thermo-mechanical  stress  analysis  can  be  performed  either  fully  coupled,  or  one-way 
coupled  in  which  only  one  disciplinary  response  affects  the  other.  Regardless  of  the 
formulation,  the  governing  equation  may  be  expressed  as 


<J  +b  =  p- 

u.i  j  r 


d2u 


(30) 


where,  <7  represents  the  stress  tensor,  b  the  body  force  terms  per  unit  mass,  p  the  mass 
density,  and  u  the  displacement  field.  To  cast  the  above  equations  in  terms  of 
displacements,  the  relationships  between  strain-displacement  and  stress-strain  must  be 
assumed.  For  small  strains  and  displacements  (i.e.,  geometric  linearity),  the  strains  are 
related  to  the  deformation  gradients  as 


£  =  —  [u  +u  ) 

v  9  V  J.‘  ‘>J ) 


(31) 


Additionally,  under  the  assumption  of  linear  elasticity  (i.e.  material  linearity),  the 
stress  may  be  related  to  strain  as 

<32> 

where,  a  is  the  coefficient  of  thermal  expansion  (which  may  in  general  be  a  function  of 
temperature),  AT  the  temperature  difference  from  reference,  and  <5  the  Kronecker  delta 
symbol.  Here,  the  first  term  relates  the  stress  to  mechanical  strain,  and  the  second  to  the 
thermal  strain.  Assuming  isotropic  material  behavior,  the  constitutive  (elasticity)  tensor 
may  be  conveniently  written  as 


Ev 

MH 


8  8k, 

ij  kl 


(33) 


where,  in  general,  the  modulus  of  elasticity,  E ,  and  Poisson’s  ratio,  v ,  may  be  functions 
of  temperature.  Currently,  it  is  assumed  that  the  mechanical  response  does  not  alter  the 
temperature  distribution,  and  therefore,  a  one-way  coupling  is  utilized.  Furthermore,  a 
steady-state  temperature  field  is  applied  to  the  structure  and,  hence,  inertia  may  be 
neglected  in  the  problem  formulation. 


The  thermo-elastic  structural  analysis  is  performed  using  a  standard  displacement- 
based  Galerkin  formulation.  With  introduction  of  the  stress-strain,  and  strain- 
displacement,  relations  into  the  equations  of  equilibrium,  the  Navier-displacement 
equations  may  be  written  (neglecting  inertia)  as 


C  ,,u,  ,.  +  b  =  0 

ijkl  kji  j 


(34) 
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Integrating  these  equations  over  the  volume  of  an  element,  using  a  standard 
Galerkin  formulation,  and  assembling  the  element  equations  yields  an  algebraic  system  to 
be  solved  for  the  unknown  nodal  displacement  vector  {d}  as 


MM=M  («) 

where,  the  global  stiffness  matrix  and  load  vector  are 

[K]  =  lKBj[c][B]dV  (36) 

M  =  XJ[>]'  W^+XjW  {'}<S+XJM’ [c]{a-AT}dV  (37) 


Here,  [iV  J  represents  the  element  shape  function  matrix,  which  relates  the 

displacement  field  to  the  element  nodal  displacements,  and  [5]  is  the  so-called  strain- 

displacement  matrix,  which  relates  the  element  strains  to  the  nodal  displacements.  The 
first  term  of  the  load  vector  represents  body  forces  acting  over  the  volume  of  the  domain, 
the  second  term  represents  forces  due  to  traction  or  stress  (/  =  0,7}.,  where  1}  are  the 

components  of  the  unit  outward  pointing  normal)  acting  on  the  boundaries  of  the  domain, 
whereas  the  last  term  are  the  loads  resulting  from  the  thermal  strains. 

The  solution  to  the  resulting  system  may  be  accomplished  with  either  an  iterative 
or  direct  method.  Iterative  solution  algorithms  can  be  beneficial  for  very  large  systems, 
which  may  prohibit  direct  solution  methods,  or  for  systems,  which  do  not  have 
symmetric,  damping  matrices.  The  direct  method  uses  a  sparse  Cholesky  decomposition 
based  on  skyline  storage  scheme.  A  preconditioned  GMRES  [64]  is  utilized  to  solve  the 
system  iteratively. 

Once  the  nodal  displacements  have  been  determined,  the  strain  and  stress  fields 
may  be  computed.  In  multidimensional  problems,  the  individual  components  of  the  stress 
tensor  do  not  provide  adequate  information  in  order  to  ascertain  the  proximity  to  failure 
or  yielding  of  the  material.  In  these  regards,  different  stress  measures  are  typically  used. 
For  brittle  materials,  the  maximum  in  plane  principal  stress  is  typically  used.  The 
principal  planes  are  those  in  which  the  shear  stress  vanishes  and,  therefore,  are  the  planes 
that  have  maximum  normal  stress.  Since  brittle  materials  tend  to  fail  due  to  normal  stress, 
appropriate  failure  criteria  are  usually  based  on  maximum  principal  stress.  For  ductile 
materials,  which  tend  to  fail  in  shear,  typically  either  the  von  Mises  or  the  Tresca  yield 
criteria  are  used.  The  von  Mises  yield  criterion  uses  the  assumption  that  the  onset  of  yield 
is  based  on  the  second  deviatoric  stress  invariant.  The  Tresca  yield  surface  is 
circumscribed  by  the  von  Mises  yield  surface,  representing  a  more  conservative  criterion 
for  prediction  of  plastic  yielding. 

3.  Design  and  Sensitivity  Analysis 

The  three-dimensional  code  developed  is  capable  of  computing  sensitivity 
derivatives  that  can  be  utilized  in  optimization  process  for  minimizing  a  specified  cost 
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function,  which  is  indicative  of  the  performance  of  the  system.  The  code  is  capable  of 
computing  sensitivity  derivatives  using  either  discrete  adjoint  method  or  direct 
differentiation  method.  Depending  upon  the  design  problem,  one  or  other  method 
provides  higher  efficiency.  For  example,  discrete  adjoint  method  is  computationally 
efficient  for  design  problems  with  large  number  of  design  variables. 

A  general  optimization  procedure  begins  by  first  defining  a  meaningful  cost 
function  and  a  desired  set  of  design  variables.  A  numerical  analysis  of  the  baseline 
system  is  then  performed.  The  results  of  the  analysis  include  the  solution  variables  Q  of 
the  discretized  partial  differential  equations,  which  are  subsequently  used  to  determine 
the  initial  cost.  Because  the  numerical  analysis  involves  discretization  of  the  partial 
differential  equations  on  a  computational  mesh,  it  should  be  noted  that  Q  represents  the 

vector  of  solution  variables  where  each  element  of  the  vector  is  representative  of  one  or 
more  physical  variables  located  at  each  mesh  point,  % . 

3.1  Discrete  Adjoint  Method 

The  cost  function  may  have  an  explicit  dependence  on  the  vector  of  design 
variables,  f. 3 ,  but  will  also  have  an  implicit  dependence  because  Q  and  X  may  also 
depend  on  the  design  variables.  Therefore,  the  cost  function  is  typically  written  to 
indicate  the  implicit  and  explicit  dependence  on  the  design  variables  as, 

f=f(Q(P),X(P),P)  (38) 

If  R  represents  the  vector  of  discrete  residuals  at  each  mesh  point,  an  augmented 
cost  function  L  can  be  defined  in  terms  of  the  original  cost  function  and  the  vector  of 
discrete  residuals  as, 

L(Q(P),X(P),P.A)=  f(Q(P),X(P),P)+  ArR(Q(P),X(P\P)  (39) 


In  equation  (39),  A  is  the  vector  of  Lagrange  multipliers  (also  known  as  costate 
variables).  Note  that  the  augmented  cost  function,  L,  is  a  scalar  quantity  that  is  identical 
to  the  original  cost  function  / ,  when  R(Q)  is  zero,  indicating  that  the  steady-state 
solution  is  obtained.  Differentiating  the  augmented  cost  function  with  respect  to  each  of 

the  design  variables  yields  the  following  set  of  equations  for-^ ,  which  is  a  column 

dp 

vector  where  each  element  represents  the  derivative  of  the  augmented  cost  function  with 
respect  to  a  particular  design  variable. 
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(40) 


Because  the  elements  of  A  are  arbitrary,  the  second  term,  which  involves  the 
derivatives  of  the  dependent  variables  with  respect  to  the  design  variables,  can  be 
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eliminated  by  solving  a  linear  system  of  equations  for  the  costate  variables,  also  known 
as  the  adjoint  equation. 


dR 

dQ 


T 


(41) 


Once  the  costate  variables  are  obtained,  the  derivatives  of  the  cost  function  with 
respect  to  all  the  design  variables  are  obtained  using  a  matrix-vector  multiplication. 
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In  numerical  simulations,  the  largest  computational  cost  of  computing  sensitivity 
derivatives  using  the  adjoint  equations  is  due  to  the  solution  of  the  analysis  and  adjoint 
equations,  both  of  which  are  independent  of  the  number  of  design  variables.  The  only 
dependency  on  the  number  of  design  variables  is  in  the  evaluation  of  equation  (42), 
which  is  generally  much  cheaper  to  compute  than  either  the  analysis  or  adjoint  solutions. 

Note  that  the  terms  in  equations  (40)  -  (42)  involve  differentiation  of  the  discrete 
residual  R  ,  the  cost  function  / ,  and  the  computational  mesh  %  with  respect  to  the 
dependent  variables  Q ,  the  design  variables  /?,  and  the  location  of  the  mesh  points^  • 
Correct  implementation  of  this  procedure  can  be  extremely  tedious  to  accomplish  by 
hand  and  the  resulting  code  can  be  difficult  to  maintain.  To  overcome  the  difficulties 
associated  with  hand  differentiation,  the  complex-variable  technique  of  Burdyshaw  et  al. 
[16]  and  Nielsen  et  al.  [24]  has  been  used  for  evaluating  all  the  terms  in  the  matrices 
required  for  solving  the  adjoint  equations  and  for  evaluating  equation  (42)  once  the 
costate  variables  have  been  obtained.  Although  the  original  development  of  the  complex- 
variable  technique  is  described  in  the  literature  [14],  a  detailed  derivation  of  complex 
variable  technique  is  also  given  by  Kapadia  et  al.  [26]  along  with  the  detailed  discussion 
on  relative  benefits  and  drawbacks  of  complex  variable  method  in  comparison  to 
automatic  differentiation  and  finite  difference  methods. 


3.2  Direct  Differentiation 

Sensitivity  derivatives  can  also  be  computed  using  the  direct  differentiation 
method.  Derivation  of  this  method  using  the  chain  rule  is  shown  in  equations  (43)  —  (46). 

dp  dp  dQ  dp  dx  dp 

Now,  R(Q(P),X(P),P)=0 

dR  dR  dR  dQ  dRdx  n 

=> — = — + - — + - —  =0 

dp  dp  dQ  dp  dxdp 


(43) 

(44) 

(45) 
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(46) 


I  BP  J  [W)  ZXW 

As  seen,  computation  of  3(7/ 3/3  is  an  essential  component  of  this  method,  which 
requires  the  solution  of  a  linear  system  of  equations  for  each  design  variable.  This 
requirement  makes  direct  differentiation  methods  computationally  expensive  for 
problems  with  many  design  variables.  However,  as  derivatives  of  dependent  variables 
with  respect  to  the  design  variables  are  computed  at  each  node  in  the  flowfield,  this 
method  is  particularly  useful  when  there  are  many  flowfield  constraints.  Direct 
differentiation  method  is  utilized  in  this  report  for  computing  sensitivity  derivatives  in  a 
case  describing  fluid-structure  interaction  capability  for  planar  SOFC. 

3.3  Mesh  Sensitivity 

Shape  is  one  of  the  most  important  parameters  in  computational  design  of 
physical  systems.  The  three-dimensional  code  is  capable  of  computing  sensitivity 
derivatives  of  desired  cost  function  with  respect  to  different  shape  parameters.  Shape 
parameter  can  range  from  a  single  mesh  point  to  the  surface  shape  and  size  of  different 
components  of  the  system.  A  parameterization  technique  to  represent  the  shape  of  the 
computational  domain  has  been  developed  and  described  in  Section  3.4.  To  maintain  the 
quality  of  the  mesh  during  a  design  cycle,  a  methodology  is  required  to  compute  the 
displacements  of  the  interior  nodes  when  the  underlying  geometry  is  modified.  The 
present  simulations  use  the  linear  elasticity  equations  as  applied  in  reference  [64]  to 
compute  these  displacements  as  shown  in  equation  (47). 

[T]X=Xsulface  (47) 

The  matrix,  [T],  is  formed  by  applying  a  finite-volume  method  to  the  linear 
elasticity  equations  and  Xsurface  denotes  the  displacements  applied  to  the  surface  nodes. 

Note  that,  [  T  ]  does  not  depend  on  the  vector  of  the  design  variables,  J3 .  Thus,  by 
differentiating  equation  (47)  with  respect  to  ft ,  mesh  sensitivities,  3// 3/?  can  be 
obtained. 

[r]d*=  ^  (48) 

L  J  3/3  3/3  v  ' 

Using  equation  (48),  mesh  sensitivities  are  computed  separately  for  each  shape 
parameter  and  are  then  used  in  equation  (42)  for  determining  the  sensitivity  derivatives  of 
the  overall  cost  function.  Because  this  procedure  is  repeated  for  each  design  variable,  it 
can  be  computationally  prohibitive  for  three-dimensional  design  problems  when  many 
parameters  are  present.  To  overcome  this  difficulty,  the  method  developed  by  Nielsen 
and  Park  [65]  has  been  implemented  in  the  current  study.  In  this  technique,  satisfaction 
of  the  mesh  equation  given  by  equation  (47)  is  included  as  a  further  constraint  in  the 
augmented  cost  function. 
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(49) 


L(Q(P),  X(P),  A  A)  =  f(Q(p),  X(P),  P)  +  A7  R(Q(P),  X(P X  £) 
+  As([r]*-*»*J 


Here,  A  is  the  vector  of  co-state  variables  associated  with  the  mesh 

g 

displacements.  The  last  term  in  equation  (49)  represents  the  residual  of  the  linear  system 
presented  in  equation  (47),  which  is  zero  when  the  solution  is  converged.  Thus,  equation 
(49)  maintains  the  original  value  of  the  desired  cost  function,  / .  Equation  (50)  is 
obtained  by  following  the  same  procedure  used  in  deriving  equation  (42). 
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(50) 


Finally,  the  grid  adjoint  problem,  equation  (51),  is  derived  by  solving  for  A  .to 

eliminate  the  mesh  sensitivity  term,  (3^7  3/?)7  in  a  similar  manner  as  the  first  adjoint 
problem. 


(51) 


Note  that  with  this  procedure,  A  is  first  obtained  and  subsequently  used  on  the 
right-hand  side  of  equation  (51).  Although  equation  (51)  represents  an  additional  linear 
system  of  equations,  the  effects  of  mesh  sensitivities  for  each  design  variable  are 
accounted  for  in  a  programming  loop  extending  only  over  surface  coordinates  and 
eliminates  the  need  for  multiple  solutions  of  equation  (47).  By  combining  equations  (41), 
(50)  and  (51),  sensitivity  derivatives  of  an  augmented  cost  function  can  be  computed 
using  equation  (52). 
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3.4  Parameterization 
3.4.1  Control  Grids 

A  parameterization  method  has  been  developed  to  improve  the  flexibility  and 
speed  in  which  a  shape  design  problem  can  be  defined.  This  method  uses  a  construct 
called  a  control  grid  [66],  which  is  associated  with  the  surface  mesh  upon  which  shape 
modification  is  desired.  Referring  to  Figure  3(a),  design  variables  are  defined  on  the 
boundaries  of  the  control  grid  as  perturbation  sources,  which  are  then  propagated  through 
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the  domain  via  the  solution  of  an  elliptic  PDE  (Laplacian)  over  the  control  grid  volume. 
The  perturbations  at  the  surface  grid  points  are  then  obtained  by  interpolation  from  the 
surrounding  nodes  on  the  control  grid.  This  procedure  is  applicable  when  disparate 
meshes  have  been  generated  for  multidisciplinary  analysis.  Because  the  resulting  surface 
displacements  are  linear  functions  of  the  design  variables,  the  parameterization  need  only 
be  computed  for  the  original  surface  mesh.  Subsequent  shape  deformations  are  then 
computed  as  a  linear  combination  of  the  design  variable  values  and  their  associated 
sensitivity  derivatives,  which  are  also  computed  only  for  original  surface.  This  tool  has 
been  used  to  define  a  parameterization  for  a  wide  variety  of  shapes,  including 
turbomachinery  blades,  wing/spar  combinations,  and  an  inlet  s-duct.  An  example  design 
case  performed  using  this  technique  is  shown  in  Figure  3(b)  and  3(c).  As  seen  in  the 
figures,  a  surface  for  an  aerodynamic  analysis  and  an  interior  mesh  representing  the 
underlying  structure  are  both  represented.  As  the  shape  changes  in  response  to 
adjustments  in  the  design  variables,  both  meshes  deform  in  unison  thereby  maintaining  a 
watertight  geometry.  This  parameterization  technique  has  been  utilized  to  provide  design 
variables  to  control  the  width  of  the  17  fuel  channels  of  a  manifold  while  treating  the 
baffles  as  rigid  bodies  [27].  Same  technique  is  applied  for  a  single-channel  SOFC,  where 
design  variables  are  located  on  the  top  and  bottom  walls  of  the  fuel  channel  [27]. 


Figure  3.  Parameterization 


3.4.2  CAPRI 

CAPRI  [67]  is  a  methodology  developed  by  the  Aerospace  Computational  Design 
Laboratory  at  the  Massachusetts  Institute  of  Technology  (MIT).  This  technology  provides 
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an  Application  Programming  Interface  (API)  to  allow  vender  neutral  communication  with 
CAD  software.  Through  coordination  with  Dr.  Bob  Haimes  at  MIT  (the  developer  of 
CAPRI),  a  framework  has  been  developed  for  integrating  CAPRI  with  the  SimCenter's 
geometry  libraries.  With  this  capability,  design  variables  used  in  the  construction  of  the 
underlying  geometry  can  be  exposed  through  the  API's,  thereby  allowing  them  to  be  used 
directly  as  design  variables  during  an  optimization  process.  The  advantage  of  this 
approach  is  that  after  a  simulation-based  design  improvement  has  been  completed;  the 
geometry'  is  represented  in  CAD  format  so  the  new  configuration  can  be  manufactured 
without  first  having  to  "re-loft"  the  geometry. 

Another  important  application  of  this  capability  lies  in  higher-order  finite  element 
simulations.  While  this  interface  had  been  originally  developed  for  purposes  of  design 
optimization,  it  is  also  used  for  placing  additional  nodes  or  quadrature  points  onto  the 
actual  surface  geometry  as  defined  by  the  CAD  definition.  During  this  process,  as 
boundary  elements  are  curved  to  conform  to  the  original  geometry  configuration, 
collapsed  elements  are  likely  to  be  generated,  particularly  when  highly  stretched  elements 
are  applied  in  the  viscous  boundary  layer.  In  this  context,  a  robust  mesh  movement 
strategy  must  be  employed  to  accommodate  the  projection  of  the  surface  meshes  for  two 
or  three-dimensional  geometries.  Here  exclusive  use  is  made  of  a  modified  linear 
elasticity  theory,  which  assumes  that  the  computational  mesh  obeys  the  isotropic  linear 
elasticity  relations.  More  information  about  this  methodology  is  given  in  Section  3.5. 


3.5  Mesh  Movement 

One  of  the  most  important  steps  in  the  shape  design  process  is  to  utilize  the 
surface  derivatives  and  modified  design  parameters  to  perform  volumetric  mesh 
deformation.  The  volumetric  shape  deformation  is  computed  by  solving  linear  elasticity 
equations,  (53)  -(55).  A  finite  volume  linear  elasticity  solver  is  utilized  to  perform  this 
task.  A  newly  developed  finite  element  method  can  also  be  used  and  supports  higher- 
order  spatial  accuracy,  which  adds  more  flexibility  and  is  compatible  with  future 
development  plans  at  the  SimCenter.  This  solver  is  utilized  in  results  shown  in  Section 
4.6. 
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In  equations  (53)  -  (55),  displacements  in  x,  y  and  z  directions  are  denoted  by  u,  v 
and  w,  respectively.  “E”  and  “v”  denote  modulus  of  elasticity  and  Poison’s  ratio, 
respectively.  While  utilizing  CAPRI  interface  and  linear  elasticity  solver  to  generate 
higher-order  meshes,  modulus  of  elasticity  is  set  to  be  proportional  to  the  inverse  of  wall 
distance,  to  enable  the  nodes  closer  to  boundaries  to  present  higher  rigidity  for 
maintaining  the  fidelity  of  the  boundary-layer  region.  The  value  of  Poison's  ratio  should 
be  in  the  range  of -1  and  0.5.  For  higher-order  finite  element  work,  a  constant  value  of 
0.25  is  set  and  is  found  to  be  effective  to  transmit  surface  deformation  into  the  interior 
while  capable  of  producing  valid  high-aspect-ratio  cells. 

3.6  Optimization 

Sensivity  analysis  can  be  used  to  choose  the  design  parameters  that  have  strong 
influence  on  objective  function  and  ignore  non-important  parameters.  This  information 
is  helpful  prior  to  an  optimization  study.  In  addition,  sensitivy  derivitives  are  essential 
part  of  many  optimization  algoithms.  especially  gradeint-based  methods.  In  gradient- 
based  methods,  gradients  of  the  response  functions  are  computed  to  find  the  direction  of 
improvement.  Gradient-based  optimization  is  the  search  method  that  underlies  many 
efficient  local  optimization  methods.  A  drawback  to  this  kind  of  methods  is  that 
computing  gradients  is  expensive.  However,  by  using  methods  like  discrete  ajoint  and 
direct  diffrention  this  cost  can  be  decreased  significantly  when  a  large  number  of  design 
parameters  are  present. 

Optimization  is  the  minimization  or  maximization  of  a  cost  function  subjected  to 
constrains  on  its  variables.  The  optimization  problem  can  be  defined  as, 

min  f(x)xemn  subject  to  gi(x)  <  0  iEl  (56) 

where  f  is  the  objective  function,  a  function  of  x  that  we  want  to  optimize,  x  is  a 
vector  of  design  variables.  There  are  several  methods  for  solving  an  optimization 
problem.  Table  3  summerizes  some  of  these  methods  and  their  properites. 

Several  optimization  software  packages  have  been  developed  in  the  past.  For 
this  study,  DAKOTA  toolkit  is  used.  The  DAKOTA  (Design  Analysis  Kit  for 
Optimization  and  Terascale  Applications)  was  developed  by  the  engineers  at  the  Sandia 
National  Laboratories  [56].  DAKOTA’S  optimization  capabilities  include  a  wide  variety 
of  gradient-based  and  nongradient-based  optimization  methods.  It  includes  many 
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extenal  optimization  libraries  such  as  the  OPT++  library  [68],  CONMIN  and  DOT 
libraries  [69]. 


Table  3.  The  optimzation  methods  and  properties 

Optimization  method 

Properties 

Gradient  descent  method 

Slow  convergence  (convergence  rate  is  linear) 

Zig-zag  back  and  forth  behavior 

Newton’s  method 

Requires  computation  of  Hessian 

Quadratic  convergence  rate 

Requires  a  good  initial  gauss 

Can  be  unstable 

Can  be  converged  to  maximum 

Gauss-Newton  method 

Can  be  unstable 

Conjugate  gradient 
method 

Requires  line  search 

Does  not  need  explicit  second  derivatives 

An  interface  is  created  to  link  the  flow  solver  to  DAKOTA.  Figure  4  shows  the 
workflow  involving  an  in-house  solver  and  DAKOTA.  The  output  file  from  DAKOTA 
containing  new  values  for  the  solver  is  "params.in”.  This  file  is  created  by  DAKOTA 
during  each  design  cycle. 


Read  initial  file  including  values  and 
constraints  for  design  variables 


Read  the  values  for  objective  function  and  ^ _  Create  the  file  output  and  write  the  cost 

gradients  from  the  file  output  function  and  gradients  to  this  file 


Figure  4.  The  data  interchange  and  interface  between  the  flow  solver  and  DAKOTA 
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3.7  Solution  Procedure 

In  the  three-dimensional  solver,  flowfield  variables  are  computed  using  an 
unstructured,  implicit,  finite-volume  scheme.  The  solver  is  vertex  centered  and  the 
discrete  residual  at  each  node  is  computed  by  integrating  the  governing  equations,  (1)  - 
(6)  over  a  median  dual  control  volume.  Because  a  steady-state  solution  is  the  primary 
goal  of  the  current  work,  time  accuracy  of  the  solution  is  sacrificed  by  allowing  local 
time-stepping  to  accelerate  convergence. 

To  reduce  computer  time,  the  solution  is  obtained  using  multiple  processors 
utilizing  the  message-passing  interface  (MPI)  [70]  and  necessary  grid  decomposition  is 
achieved  using  METIS  [71].  Original  grids  are  generated  using  the  commercial  software 
Pointwise  [72]. 

An  implicit  Euler  scheme  is  used  to  solve  the  non-linear  system  as  given  by 
equations  (1)  -  (6).  A  flux-difference  splitting  scheme  based  on  the  ROE  scheme  [73,74] 
for  a  multi-component  mixture  is  derived  to  model  the  convective  fluxes.  A  central- 
difference  formulation  is  used  to  compute  all  the  second-order  derivative  terms.  Linear 
systems  encountered  in  both  the  flowfield  and  sensitivity  solvers  are  solved  using  the 
GMRES  [64]  method  with  ILU-K  preconditioner. 

A  domain  decomposition  tool  has  been  developed  to  split  computational  domain 
into  any  number  of  subdomains.  This  tool  is  capable  of  generating  internal  interfaces, 
which  are  essential  for  multidisciplinary  applications  that  contain  domains  with  different 
properties.  The  tool  is  also  capable  of  generating  higher  order  grids  using  initial  linear 
grid  for  higher-order  finite  element  applications.  Currently,  this  tool  is  being  utilized  for 
multidisciplinary  simulation  applications  including,  CFD,  electromagnetics  and  batteries. 

4.  Results  and  Discussion 

4.1  Code  Validation 

In  this  section,  code  validation  results  for  baseline  solution  and  sensitivity 
derivatives  are  presented  for  SOFC  simulations.  The  geometry  utilized  in  these 
simulations  is  the  same  as  used  by  Wang  et  al.  [57].  One  of  the  main  reasons  for  using 
this  geometry  is  the  availability  of  an  experimental  polarization  curve  [57]  that  can  be 
used  to  validate  the  in-house,  three-dimensional  multi-species  Navier-Stokes  solver. 

4.1.1  Baseline  Solution 

The  fuel  mixture  is  assumed  to  contain  hydrogen  and  steam  only,  i.e.  chemistry  is 
not  present  in  the  validation  case.  Air  is  modeled  as  a  mixture  of  oxygen  and  nitrogen. 
Species  mole  fractions  of  the  fuel  mixture  and  air  entering  the  respective  channels  are 
given  in  Table  4.  The  operating  pressure,  temperature  and  mass  flow  rates  of  fuel  and  air 
are  also  given  in  Table  4.  Wang  et  al.  [57]  has  obtained  the  polarization  curve  under  the 
same  operating  conditions  as  described  in  Table  4.  Surface  grid  of  the  computational 
domain  is  shown  in  Figure  5. 
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Table  4.  Operating  conditions  utilized  in  polarization  curve 

** 

** 

T(K) 

P(N!m2) 

g 

m  fuel  (kg  /  s) 

g 

ma,r  (kg/s) 

0.9578 

0.0422 

0.198 

0.802 

1273  K 

101325 

5.94X10-7 

2.15X10-5 

A  comparison  between  the  experimental  polarization  curve  [57]  and  the 
polarization  curve  obtained  using  the  numerical  model  is  shown  in  figure  6.  As  can  be 
seen,  the  overall  comparison  is  satisfactory.  The  numerical  tool  successfully  predicts  the 
shape  of  the  polarization  curve  and  obtains  results  that  are  within  two  percent  of  the 
experimental  data  at  low  current  densities  and  are  essentially  identical  to  the  experimental 
data  at  higher  current  densities.  As  expected,  the  cell  voltage  reduces  with  increasing 
current  density  due  in  part  to  Ohmic  losses  which  are  linearly  proportional  to  the  current 
density.  Also,  increase  in  current  draws  more  hydrogen  and  oxygen  from  the  anode- 
electrolyte  and  the  cathode-electrolyte  interfaces,  respectively  and  produces  more  steam. 
This  reduces  the  value  of  “EMF”  in  equation  (18)  i.e.  concentration  polarization 
increases.  Thus,  the  cumulative  effects  of  Ohmic  polarization  and  concentration 
polarization  are  evident  in  Figure  6. 


Figure  5.  Surface  Mesh 
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Figure  6.  Polarization  curve 


4.1.2  Sensitivity  Derivatives 

In  this  study,  the  following  two  cost  functions  are  considered. 

Cost  -1 :  Cell  voltage  —  equation  (57.1). 

Cost  -2:  Term  responsible  for  the  concentration  polarization  at  the  anode-electrolyte 
interface  -  equation  (57.2). 


It  should  be  noted  that  the  sensitivity  results  described  in  this  section  assume 
H2  /  H20  as  the  fuel  mixture.  The  operating  conditions  are  the  same  as  described  in 

Table  4  and  the  applied  current  density  is  4000  A I  m1 . 

Improving  power  output  is  the  ultimate  goal  of  the  SOFC  design.  If  current 
density  is  fixed,  the  power  output  can  be  improved  by  increasing  the  cell  voltage.  Thus, 
the  first  cost  function  chosen  in  this  study  is  the  cell  voltage.  Sensitivity  derivatives  of  the 
cost  function  representing  the  cell  voltage  with  respect  to  various  design  parameters  can 
be  extremely  useful  in  the  design  cycle. 
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The  second  cost  function  represents  the  term  responsible  for  the  concentration 
polarization  at  the  anode-electrolyte  interface.  As  seen,  equation  (57.2)  exhibits  an 
interesting  nature  due  to  its  explicit  dependence  on  the  species  partial  pressures  and  the 
temperature.  Also,  improvement  in  concentration  polarization  can  increase  the  SOFC 
performance  and  thus,  it  is  chosen  as  the  second  cost  function  in  sensitivity  studies. 

Six  design  parameters  are  included  to  compute  sensitivity  derivatives  of  the 
aforementioned  objective  functions.  The  design  parameters  are  comprised  of  the  material 
properties  of  the  anode  and  the  cathode  including  porosity,  tortuosity  and  mean  pore 
radius. 


To  validate  the  implementation  of  the  discrete  adjoint  method,  sensitivity 
derivatives  computed  using  the  discrete  adjoint  method  are  compared  with  the  same 
derivatives  computed  using  the  direct  differentiation  method  and  the  finite  difference 
method.  Note  that  the  sensitivity  derivatives  computed  using  central  finite  difference 
method  require  two  flowfield  solutions  for  each  design  variables.  To  reduce 
computational  efforts,  the  comparison  study  is  performed  using  a  single  channel 
geometry  (original  geometry  contains  six  channels)  and  coarse  grid.  Relevant  physics 
included  for  the  original  geometry  is  included  in  the  comparison  study. 

Tables  5  and  6  show  the  comparisons  amongst  sensitivity  derivatives  computed 
using  the  different  methods  for  the  Cost-1  and  Cost-2,  respectively.  As  seen,  sensitivity 
derivatives  computed  using  the  discrete  adjoint  method  and  direct  differentiation  method 
match  up  to  9  to  11  digits.  This  comparison  is  excellent.  Also,  matching  significant 
digits  between  finite  difference  results  and  the  discrete  adjoint  method  results  vary 
between  2  to  4.  Due  to  subtractive  cancellation  errors,  it  is  hard  to  find  an  optimum  step 
size  when  the  finite  difference  method  is  used  to  compute  sensitivity  derivatives.  Thus, 
comparison  between  finite  difference  derivatives  and  the  discrete  adjoint  derivatives  are 
considered  satisfactory. 


Table  5.  Validation  of  adjoint  implementation  (Cost  1) 

D.V. 

Discrete  Adjoint 

Direct  Differentiation 

Finite  Difference 

£a 

-1.4136629036e-02 

-1.4136629036e-02 

-1.411798953e-02 

Ka 

-3. 4 92 4988954 e- 03 

-3. 4 92 4988954 e- 03 

-3. 49100096 6 e-03 

<r>a 

8. 757810 687 0e+02 

8. 7578106869e+02 

8.754119812e+02 

£c 

2.7292696323e-03 

2 . 72 92 69 632 2e- 03 

2 . 761454211e-03 

K 

-1 . 4 9760417 63e-03 

-1 . 4 97 60 4 17 63e- 03 

-1. 4973217 55e-03 

<r>c 

1.8945315028e+02 

1. 8945315028e+02 

1.894163071e+02 

A  run-time  comparison  between  the  discrete  adjoint  and  direct  differentiation 
method  is  shown  in  Table  7.  Note  that  this  comparison  is  made  using  the  coarse  mesh 
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utilized  for  Tables  5  and  6.  As  mentioned  before,  the  direct  differentiation  method 
requires  the  solution  of  a  linear  system  for  each  design  variable.  Thus,  the  run-time  of  a 
direct  differentiation  method  is  a  linear  function  of  the  number  of  design  variables.  On 
the  other  hand,  the  discrete  adjoint  method  only  requires  the  solution  of  a  single  linear 
system  and  the  dependency  on  the  number  of  design  variables  is  that  a  matrix- vector 
product  must  be  computed  for  each  design  variable. 


Table  6.  Validation  of  adjoint  implementation  (Cost  2) 

D.V. 

Discrete  Adjoint 

Direct  Differentiation 

Finite  Difference 

£a 

-1. 387075190 6e-02 

-1.3870751907e-02 

-1 . 386688  67e-02 

K 

-3 . 5 95 9382 4 5 Oe- 03 

-3 . 5 95 93 82 4 50e- 03 

-3.59483888e-03 

<r>a 

9 . 0264297426e+02 

9 . 0264297426e+02 

9 . 0237  9501e+02 

1. 6938729757 e- 03 

1. 6938729757e-03 

1 . 6919354 le-03 

Kc 

8.7332085595e-05 

8 . 7 332 08 55 95e- 05 

8 . 754  93837e-05 

<r>c 

-1 . 01 643125 15e+01 

-1.0164312515e+01 

-1 . 01922398e+01 

Table  7.  Run-time  comparison  between  adjoint  and  direct  differentiation  methods 

Adjoint 

Direct  Differentiation 

Run-time  (sec) 

125 

571 

Once  implementation  of  the  discrete  adjoint  method  is  validated,  the  next  step  in 
the  design  process  is  to  compute  sensitivity  derivatives  for  the  actual  geometry.  Table  8 
shows  the  sensitivity  derivatives  of  both  cost  functions  obtained  using  the  discrete  adjoint 
method  for  the  original  geometry. 


Table  8.  Sensitivity  derivatives  computed  using  discrete  adjoint  method 

D.V. 

A 

A 

-1. 0110037400 e- 02 

-1. 0899218104 e- 02 

-5.0238174308e-03 

-5.1098651595e-03 

<r>a 

1. 1323025342 e+ 03 

1 . 154  94  9837  0e+03 

3.3425032057e-03 

2.3809517707e-04 

-1. 6269401390e-03 

-1.29598 07 45 Oe- 05 

<r>c 

2. 0610982993 e+02 

-7 .5667688361e-01 
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4.2  Mass  Flow  Rate  Effects 

In  this  section,  effects  of  mass  flow  rate  variations  on  objective  function  are 
investigated  for  two  different  cases.  The  first  case  utilizes  pure  hydrogen  as  fuel  and  the 
second  case  utilizes  fuel  mixture  containing  methane.  The  second  case  also  includes 
chemical  reactions  of  steam  reforming  as  well  as  shift  reaction. 


Mass  flow  enforcement  is  required  to  accurately  model  the  physics  of  the  given 
problem.  Previously,  continuous  back  pressure  adjustment  was  required  to  achieve  the 
desired  mass  flow  rate.  Recently,  better  control  over  mass  flow  rate  enforcement  has  been 
achieved  by  re-formulating  boundary  conditions  including,  inlet,  outlet  and  no-slip  walls. 


Figure  7  shows  the  relationship  between  the  mass  flow  rate  of  fuel  and  the  cell 
voltage.  The  fuel  composition  and  operating  conditions  (temperature  and  pressure)  are 
the  same  as  shown  in  Table  4.  The  air  mass  flow  rate  is  also  constant  at 


mwr  =  1 .76 x  1 0-4 kg / sec  and  the  applied  current  density  is  specified  as 4000  A/m2 .  Fuel 

mass  flow  rate  has  been  increased  gradually  to  analyze  its  effect  on  output  voltage.  As 
can  be  seen,  the  cell  voltage  increases  with  increases  in  fuel  mass  flow  rate.  As  described 
earlier,  the  three  major  contributors  to  the  voltage  loss  in  SOFC  are  Ohmic  polarization, 
concentration  polarization  and  activation  polarization.  For  the  particular  case  shown  in 
Figure  7,  concentration  polarization  plays  a  major  role  in  deciding  the  variation  in  cell 
voltage  with  the  fuel  mass  flow  rate.  Concentration  polarization  strongly  depends  on  the 
partial  pressures  of  hydrogen,  steam  and  oxygen,  as  given  by  equation  (18).  The  term 


responsible  for  the  concentration  polarization  in  equation  (18)  is  In 


f  P„  Po 5  A 

ti  9  CA 


[  H,0 


The  value 


of/g5 


can  be  assumed  constant  for  different  cases  in  Figure  7  as  current  density  and  air 

mass  flow  rate  remain  constant.  Thus,  the  determining  factor  for  the  voltage  loss  is  the 
ratio  of  partial  pressures  of  hydrogen  to  steam,  which  becomes  smaller  as  the  mass  flow 
rate  of  fuel  is  reduced  and  thus,  affecting  concentration  polarization  adversely.  As  seen 
in  Figure  7,  the  effect  of  concentration  polarization  is  weak  for  the  cases  with  higher 
mass  flow  rate  but  becomes  more  prominent  once  the  fuel  mass  flow  rate  reduces  below 

=  1 .0  x  1  (T6  kg  /  sec ,  which  indicates  that  some  part  (near  the  fuel  channel  outlet)  of 


m 


fuel 


the  anode-electrolyte  interface  may  be  starving  due  to  high  fuel  utilization.  This 
eventually  contributes  towards  a  reduction  in  electrochemical  activity  in  the  affected 
regions. 
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Figure  7.  Voltage  vs.  Mass  flow  rate  (Operating  conditions:  Table  4) 

Results  presented  in  Figure  7  utilize  a  H2  /  H20  mixture  as  fuel  and  does  not 

include  chemistry  inside  the  electrodes.  To  demonstrate  the  combined  effect  of 
chemistry/electrochemistry  on  SOFC  performance,  a  fuel  mixture  containing 
H2,H20, CO,  C02  and CH4  is  used  in  the  next  case.  Air  is  modeled  as  a  mixture  of 

02!N2.  Mole  fractions  of  various  species  and  operating  conditions  utilized  in  this  case 

are  described  in  Table  7.  Two  chemical  reactions,  namely,  methane  reforming,  equation 
(57.1)  and  water-gas  shift  reaction,  equation  (58.2)  are  present  inside  the  anode. 

CHa  +  H20  ^=>  CO  +  3 H2 

CO  +  H.O  t  —  ^  CO,  +  H? 

2  2  2 

Methane  reforming  is  an  endothermic  reaction  and  tends  to  reduce  the 
temperature  of  the  system.  The  shift  and  electrochemical  reactions  are  both  exothermic 
and  tend  to  increase  the  temperature  of  the  system.  Also,  methane  reforming  is  a  fast 
reaction  that  finishes  rapidly  in  the  region  located  near  the  fuel  channel  inlet.  The  size  of 
this  region  depends  on  the  mass  flow  rate  of  the  fuel. 


(58.1) 

(58.2) 
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Table  7.  Operating  conditions  utilized  in  Figure  8 

^co 

Xh-o 

P(N/m2) 

T(K) 

** 

0.029 

0.493 

0.044 

0.263 

0.171 

101325 

1273  K 

0.198 

0.802 

To  investigate  the  effect  of  fuel  mass  flow  rate  on  the  cell  voltage,  several  cases 
were  run  using  operating  conditions  described  in  Table  7  with  a  fixed  air  flow  rate  of 

mair  =  1 .76  x  1 0"4  kg  /  sec  and  constant  current  density  of  4000  A/m2.  The  fuel  mass  flow 

rate  is  reduced  gradually  and  its  effect  on  voltage  output  is  plotted  in  figure  8.  Even 
though  the  operating  conditions  utilized  in  the  present  case  (figure  8)  and  the  previous 
case  (figure  7,  H2  /  H20  mixture)  are  the  same  except  for  the  fuel  type,  the  nature  of  plots 
is  opposite  in  figures  7  and  8. 
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Figure  8.  Voltage  vs.  Mass  flow  rate  (Operating  conditions:  Table  7) 

As  seen  in  figure  8,  cell  voltage  reduces  with  increases  in  the  fuel  mass  flow  rate, 
which  is  opposite  to  the  trend  found  in  figure  7.  As  described  earlier,  the  reduction  in 
flow  rate  of  fuel  increases  concentration  polarization,  which  is  responsible  for  the  voltage 
drop.  In  the  previous  case,  as  no  chemical  reactions  were  present  (thus,  no  hydrogen 
source),  the  relationship  between  the  fuel  flow  rate  and  the  cell  voltage  is  straightforward. 
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However,  the  same  is  not  true  for  the  current  case  and  the  reason  is  described  in  the  next 
paragraph. 


Concentration  of  H^O  increases  at  the  anode-electrolyte  interface  due  to  the 


electrochemical  reaction  and  additional  H20  diffuses  inside  the  anode.  This  additional 
H20  triggers  methane  reforming  and  shift  reactions  described  by  equations  (58.1)  and 
(58.2),  respectively.  In  these  reactions,  H20  participates  as  reactant  and  produces 


hydrogen  and  thus,  increases  the  value  of  the  term  In 


1  H,0 


Due  to  these 


complicated  interactions  amongst  chemical-electrochemical  reactions,  the  concentration 
polarization  does  not  necessarily  increase  with  the  reduction  in  fuel  flow  rate  for  the  plot 
shown  in  figure  8.  Thus,  other  modes  of  voltage  losses  or  polarizations  are  required  to  be 
investigated  to  determine  the  reasons  behind  the  nature  of  the  plot  shown  in  figure  8. 


1273.0  1282.3  1291.7  1301.0  1310.3 

Figure  9.  Temperature  contours  at  the  anode-electrolyte  interface 
( mfuel  =  4.98x1 0~*kg/  sec) 

Figures  9  and  10  show  temperature  contours  plotted  at  the  anode-electrolyte 
interface  for  two  different  points  in  figure  8.  Figure  9  represents  the  location  with  the 

lowest  mass  flow  rate  of  fuel  (thM  =  4.98  X  lO-6^  /  sec)  and  figure  10  represents  the 
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location  with  the  highest  fuel  flow  rate (miuel  =30.02x10”*’% /sec)  .  The  difference 

between  the  lowest  and  the  highest  temperature  points  in  both  figures  is  approximately  40 
K.  However,  the  average  temperature  in  figure  9  is  approximately  73  K  higher  than  that 
in  figure  10.  Even  though  the  electrochemical  reaction  produces  the  same  amount  of  heat 
in  both  cases  (due  to  the  same  current  density),  fluid  that  is  flowing  slower  (figure  9)  is 
heated  more  than  the  fluid  that  is  flowing  faster  (figure  1 0).  As  described  earlier, 
methane  reforming  is  an  extremely  fast  endothermic  reaction  and  thus,  the  reduction  in 
temperature  caused  by  the  reforming  reaction  remains  similar  in  both  cases.  The  next 
step  is  to  determine  the  effect  of  temperature  on  Ohmic  losses  in  both  cases. 


1200  0  12094  1218  8  1228.2  1237  6 


Figure  1 0.  Temperature  contours  at  the  anode-electrolyte  interface 
{mfuel  =30.02xl0_6%/sec) 

Figures  11  and  12  show  the  electric  potential  distribution  inside  the  electrolyte  (z- 
plane)  for  the  cases  representing  the  lowest  {nifuei  =  4.98  x  10”/’%  /  sec)  and  the  highest 

( rnfud  —  30.02  xlO^&g/  sec)  mass  flow  rates  of  fuel,  respectively.  The  electrolyte 

thickness  is  0.05  mm  and  acts  as  a  main  contributor  in  Ohmic  losses.  As  seen,  the 
voltage  drop  inside  the  electrolyte  for  the  case  with  the  highest  fuel  flow  rate  (figure  1 2) 
is  almost  twice  the  voltage  drop  for  the  case  with  the  lowest  fuel  flow  rate  (figure  11). 
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Electric  Potential  (V) 


0  891  0  894  0.897  0.9  0  903  0.906 


Figure  1 1 .  Ohmic  losses  inside  the  electrolyte  (m/ue/  =  4.98 X  / sec) 


Electric  Potential  (V) 


0  868  0.877333  0  886667  0  896 


Figure  12.  Ohmic  losses  inside  the  electrolyte  {tnfuel  =30.02x10 ^kg  /  sec) 

As  mentioned  earlier,  the  electric/ionic  resistivities  of  various  components  of  the 
SOFC  are  strong  functions  of  the  temperature.  To  further  investigate  the  effect  of 
temperature,  various  values  of  electric/ionic  resistivities  for  the  anode,  the  cathode  and 
the  electrolyte  are  shown  in  Table  8  for  a  temperature  range  of  1 100  K  -  1300  K.  As 
seen,  the  ionic  resistivity  of  the  electrolyte  reduces  by  almost  half  when  the  temperature 
is  raised  from  1 100  K  to  1200  K.  The  same  is  true  when  the  temperature  is  increased 
from  1200  K  to  1300  K.  This  dependence  can  be  correlated  with  the  voltage  drop  found 
in  figures  1 1  and  12  for  the  cases  with  the  lowest  fuel  flow  rate  and  the  highest  fuel  flow 
rate,  respectively.  Note  that  the  electric  resistivity  of  the  cathode  is  also  reduced  with 
increase  in  temperature  but  that  the  increase  is  not  as  significant  as  found  for  the 
electrolyte  and  also,  the  value  of  the  resistivity  is  approximately  three  orders  smaller  than 
that  in  the  electrolyte.  Thus,  the  voltage  drop  due  to  the  Ohmic  polarization  in  the 
cathode  can  be  assumed  to  be  constant  for  both  cases.  Also,  Ohmic  polarization  inside 
the  anode  can  be  ignored  due  to  a  very  small  size  (compared  to  the  electrolyte  resistivity) 
of  the  electric  resistivity  of  the  anode. 


Table  8.  Relationship  between  electric/ionic  resistivity  and  temperature 

Temperature  (K) 

Anode 

Cathode 

Electrolyte 

1100 

8.40691e-06 

0.000139929 

0.358644 

1200 

9.34189e-06 

0.000133711 

0.163733 

1300 

1.02138e-05 

0.000128666 

0.0843334 
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4.3  Catalytic  Partial  Oxidation  of  Methane 
4.3.1  Baseline  Solution  and  Validation 

In  this  section,  catalytic  partial  oxidation  of  methane  over  RI1/AI2O3  coated 
honeycomb  is  numerically  investigated.  The  honeycomb-structured  reactors  are  widely 
used  in  many  engineering  applications  such  as  the  fuel  reformers,  catalytic  convertor 
and  gas  turbine  combustors.  For  code  validation  purposes,  experimental  work 
conducted  by  Hettel  et  al.  [50]  is  selected.  The  reactor  is  a  2  cm  diameter  cylinder  with 
260  channels  and  a  channel  density  of  600  cpsi  (channels  per  square  inch).  The  initial 
and  boundary  conditions  are  summarized  in  Table  9.  The  simulations  are  performed 
with  the  detailed  heterogeneous  oxidation  mechanism  proposed  by  Deutschmann  et  al. 
[75].  It  includes  38  heterogeneous  reactions  and  20  surface-adsorbed  species.  The  site 
density  is  assumed  to  be  2.79xl0-9  mol/cm2.  The  kinetic  data  of  the  surface-reaction 
mechanisms  are  taken  from  the  literature.  Eight  gas-phase  species  (CH4,  CO2,  H2O,  N2, 
O2,  CO,  OH  and  H2)  are  considered  for  the  simulation.  The  surface  chemistry  is 
modeled  using  the  mean-field  approximation.  Homogenous  combustion  in  the  gas  phase 
is  ignored  in  this  study  since  it  has  no  the  significant  effect  on  the  flow  field  for  this  test 
case  and  operating  conditions  [76].  Computational  grid  is  comprised  of  122208 
tetrahedral  cells.  Figure  13  shows  the  surface  grid  of  the  computational  domain 
representing  one  channel  of  the  monolith.  The  grid  is  refined  in  the  regions  near 
catalytic  wall  to  accurately  resolve  boundary  layers.  The  “inflow’'  boundary'  condition  is 
used  at  the  channel  inlet  and  a  fully  developed  boundary’  condition  is  considered  for  the 
outlet.  The  no-slip  boundary  condition  with  catalytic  reaction  source  term  is  applied  at 
the  channel  walls.  The  temperature  of  the  catalytic  wall  is  assumed  to  be  constant  along 
the  channel.  The  nonlinear  system  of  equations  obtained  from  the  discretization  is 
solved  using  Newton’s  method.  The  convergence  history'  of  the  solution  is  shown  in 
figure  14. 


Table  9.  Initial  and  boundary  conditions  for  catalytic  combustion  of  methane 

Gas  inlet  velocity 

0.329  m/s 

Gas  inlet  temperature 

1000  K 

Wall  temperature 

1000  K 

Gas  inlet  compositions(mole  fraction) 

xCHt  =  0.133 

x0  =  0.067 

N 

II 

o 

bo 

Working  pressure 

1  atm 

Channel  width 

1  mm 

Channel  length 

10  mm 
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Figure  13.  Surface  grid 


Figure  14.  Convergence  history 


Figure  15  shows  the  comparison  between  the  numerical  results  and  experimental 
data  for  the  species  mole  fractions  as  a  function  of  position  in  the  reactor.  As  indicated 
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in  the  figure  15,  the  good  agreement  is  obtained.  The  methane  oxidation  starts  instantly 
at  the  inlet  of  the  reactor.  As  expected,  oxygen  is  completely  consumed  in  the  first  few 
millimeters  of  the  reactor. 


Figure  15.  Comparison  between  the  numerical  results  and  experimental  data  for  the  partial 

oxidation  of  methane 
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Figure  16.  Species  mole  fraction  along  reformer  symmetry  axis  for  Rh  and  Pt  catalyst 
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Rhodium  and  platinum  are  considered  good  catalysts  in  terms  of  stability  and 
yields.  They  are  widely  used  for  partial  oxidation  and  catalytic  combustion  of  methane 
in  fuel  reformers,  catalytic  burners  and  catalytic  gas  turbines.  To  better  understand  the 
performance  of  methane  reformer  with  these  two  catalysts,  numerical  simulations  are 
performed.  The  detailed  heterogeneous  oxidation  mechanisms  developed  by 
Deutschmann  et  al.  [75]  (24  heterogeneous  reactions  and  1 1  surface-adsorbed  species) 
and  Deutschmann  et  al.  [77]  (38  heterogeneous  reactions  and  20  surface-adsorbed 
species)  are  used  to  model  surface  chemistry  for  rhodium  and  platinum,  respectively. 
The  temperature  of  catalyst  wall  is  fixed  to  1070  K.  The  inlet  velocity  is  considered  to 
be  0.5  m/s.  Figure  16  shows  the  mole  fraction  of  species  along  symmetry  axis  of  the 
reformer  for  both  catalysts.  As  seen,  oxygen  is  completely  consumed  (conversion  of 
99%)  in  both  cases.  Rhodium  shows  better  performance  for  partial  oxidation  of  methane 
(conversion  of  90%)  than  platinum  (conversion  of  77%). 
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(e) 


(f) 


Figure  17.  Species  mole  fractions  and  streamwise  velocity  (Pt  catalyst) 


Figures  17  and  18  show  species  mole  fraction  contours  for  reactors  with  platinum 
and  fhodium.  respectively.  Streamwise  velocity'  contours  are  also  shown  in  figure  17(f). 
The  gradient  of  hydrogen  mole  fraction  is  smaller  across  the  cross  section  of  the  channel 
as  hydrogen  has  higher  diffusion  coefficient  relative  to  other  species  considered  in  this 
simulation.  The  maximum  velocity  in  the  channel  is  close  to  1  m/s. 


<  H2 


(a) 
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1  o? 


(e) 

Figure  18.  Species  mole  fractions  (Rh  catalyst) 


4.3.2  Parameter  study 

In  this  section,  effect  of  the  different  design  parameters  on  the  fuel  reformer 
performance  is  investigated.  There  are  many  the  design  parameters  affecting  the  reactor 
efficiency  for  fuel  reforming.  These  design  parameters  can  be  related  to  the  shape/size 
of  the  reformer  as  well  as  operating  conditions  and  catalyst  material.  In  this  work,  inlet 
methane/oxygen  ratio,  inlet  velocity  and  catalytic  wall  temperature  are  considered  as 
design  variables.  By  using  inlet  velocity  as  one  of  the  parameters,  effect  of  different 
Reynolds  numbers  is  also  studied  indirectly  on  reformer  performance. 
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- —  CH4 


Figure  19.  Comparison  of  the  species  mole  fractions  along  reformer  symmetry  axis  for 

inlet  velocities  0.5m/s  and  5  m/s 

The  baseline  conditions  for  this  study  are  shown  in  Table  10.  Figure  19  shows 
the  comparison  of  the  mole  fraction  of  species  along  the  symmetry  axis  of  the  reformer 
with  two  different  inlet  velocities  of  0.5  and  5  m/s.  The  conversion  of  methane  is 
predicted  to  decrease  with  increase  in  inlet  velocity.  The  rate  of  oxygen  consumption 
along  the  reactor  is  also  decreased  and  therefore  the  peak  of  H2O  concentration  is 
shifted  towards  the  middle  of  the  channel  for  the  case  with  higher  inlet  velocity.  Mole 
fraction  contours  of  different  species  for  the  reactor  with  inlet  velocity  of  5  m/s  are 
presented  in  figure  20. 


Table  10.  Baseline  conditions  for  catalytic  combustion  of  methane 

Gas  inlet  velocity 

0.5  m/s 

Gas  inlet  temperature 

1070  K 

Wall  temperature 

1070  K 

Gas  inlet  compositions  (mole  fraction) 

*ch4  =  0-133 

x0  =  0.067 

CO 

o 

tl 

1  N 
£ 

* 

Pressure 

1  atm 

Channel  width 

1  mm 

Channel  length 

10  mm 

Catalyst 

Rh 
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(d) 


(e) 


Figure  20.  Species  mole  fractions  (V  =  5  m/s) 
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The  influence  of  catalytic  wall  temperature  variation  on  species  conversion  rates 
is  shown  in  figure  21 .  The  numerical  result  predicted  the  conversion  of  methane 
increases  from  90%  at  1070  K  to  96%  at  1170  K.  The  hydrogen  production  is  also 
increased  by  about  10%  at  higher  temperature. 


— - CH4 

- H2 


Figure  21.  Comparison  of  the  species  mole  fractions  along  reformer  symmetry'  axis  with 

different  catalytic  wall  temperature 

The  inlet  methane-oxygen  ratio  is  another  important  design  parameter.  Two 
different  cases  with  methane-oxygen  ratios  of  1  (xCH4  =  0.1  and  x02  =  0.1)  and  1/3 
(xCH4  —  0.05  and  x02  =  0.15)  are  numerically  investigated.  Results  obtained  are 
shown  in  figure  22,  which  compares  the  aforementioned  cases  with  a  baseline  case  of  of 
xCH4  =  0.133  and  x02  =  0.067.  Figure  22  shows  the  influence  of  variation  in  methane- 
oxygen  ratio  on  reformer  performance.  As  shown  in  figure  22,  size  of  active  methane 
conversion  region  increases  with  increasing  methane-oxygen  ratio  at  the  inlet.  The 
hydrogen  production  reaches  the  highest  value  for  richer  mixture  (ratio  of  2.0,  baseline 
case).  Hydrogen  production  for  the  mixture  ratio  of  1/3  is  minimal  and  CO,  CO2  and 
H2O  are  the  main  productions.  The  mole  fraction  contours  for  the  reformer  with  the 
methane-oxygen  ratio  of  1/3  are  shown  in  figure  23. 
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—  —  CH4 

—  —  H2 


Figure  22.  Effect  of  methane-oxygen  ratio  on  reformer  performance 
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Figure  23.  Species  mole  fractions  (methane-oxygen  ratio  is  1:3) 


The  type  of  catalyst  used  significantly  affects  performance  of  the  reformer.  The 
reactor  performance  for  two  different  catalysts,  platinum  and  rhodium,  was  investigated 
earlier.  In  addition  to  the  type  of  catalyst,  catalyst  loading  is  also  an  import  factor  in 
design  and  optimization  of  catalytic  reactors.  For  considering  the  effect  of  catalyst 
loading,  two  parameters  Fcat/geo  (ratio  of  catalytic  surface  area  to  geometric  surface 
area)  and  r]  (effectiveness  factor)  are  considered.  As  mentioned  earlier,  the 
heterogeneous  flux  with  catalyst  loading  effects  can  be  written  as, 


FluXf iet  —  F.at/g^riMWiSi 
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Figure  24.  Effect  of  variation  in  catalyst  loading  on  reformer  performance 

The  effect  of  the  catalyst  loading  factor  ( Fcatjgeor\ )  on  reformer  performance  is 
numerically  investigated.  As  shown  in  figure  24,  methane  conversion  increases  at  high 
catalyst  loading.  Oxygen  is  almost  completely  consumed  by  surface  reactions  in  the 
first  two-millimeter  of  the  reactor  length  for  all  three  cases.  The  methane  conversion 
increases  from  57  %  at  Fcat/geoT)  =  0.5  to  98  %  at  Fcat/geor\  =  2.  Hydrogen  production 
at  Fcat/geori  =  2  increases  about  7%  relative  to  the  baseline  case.  Species  mole 
fractions  for  the  Fcat/geoi)  =  0.5  case  are  plotted  in  figure  25. 


(a) 


(b) 
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(c) 


(d) 


(e) 


Figure  25.  Species  mole  fractions  ( Fcat/geoV  =  0-5) 

4.3.3  Sensitivity  Analysis 

In  the  previous  section,  relationship  between  reformer  performance  and  a  design 
variable  was  obtained  by  running  a  baseline  simulation,  changing  the  parameter  value 
and  re-evaluating  the  change  in  the  cost  function.  More  sophisticated  way  for  such 
analysis  and  computational  design  is  using  sensitivity  analysis.  In  this  method,  the 
sensitivities  are  obtained  by  computing  gradients  or  derivatives  of  the  objective  function 
with  respect  to  design  parameters  of  interest.  There  are  several  methods  for  computing 
sensitivity  derivatives.  In  this  work,  the  direction  differentiation  and  adjoint  methods 
are  used  to  obtain  the  sensitivity  derivatives.  Implementation  details  of  both  these 
methods  are  given  in  earlier  sections  of  this  report. 

The  mean  value  of  hydrogen  concentration  at  the  outlet  boundary  is  considered 
as  the  cost  function  as  increasing  the  value  of  this  cost  function  improves  the 
performance  of  the  reformer.  The  design  variables  considered  in  this  study  include  inlet 
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velocity,  methane  density,  oxygen  density,  catalytic  wall  temperature  and  catalytic  area 
ratio.  For  validation  purposes,  sensitivity  derivatives  are  obtained  using  both  adjoint  and 
direct  differentiation  methods  and  provided  in  Table  1 1 .  The  baseline  conditions  for  this 
study  are  shown  in  Table  9.  As  ilustrated  in  the  Table  1 1,  the  sensitivity  derivatives 
obtained  by  direct  differentiation  and  adjoint  methods  show  good  agreement. 


Table  1 1.  Comparison  of  sensitivity  derivatives  computed  using  different  methods 

Design  variable 

Direct  Differentiation 

Adjoint 

Inlet  velocity 

-0.00580817201 

-0.005840817202 

Inlet  methane  density 

0.10718769820 

0.107187698123 

Inlet  oxygen  density 

-0.07966058834 

-0.07966058831 

Catalytic  wall 
temperature 

0.00105843297 

0.00105843127 

Catalytic  area  ratio 

0.00113610487 

0.00113610454 

Sensitivity  derivatives  are  very  useful  for  design  and  optimization  purposes. 
With  these  derivatives,  one  can  estimate  how  important  a  design  parameter  is  for  the 
given  cost  function.  Thus,  parameters  that  are  not  important  can  be  ignored  during 
optimization  while  concentrating  on  important  variables.  For  example,  in  Table  11,  inlet 
methane  concentration  is  an  important  variable  for  the  cost  function  of  outlet  hydrogen 
concentration. 

4.3.4  Optimization 

As  mentioned  earlier,  an  interface  is  created  to  link  the  flow  solver  code  to 
DAKOTA  to  perform  optimization.  In  the  first  case,  numerical  optimization  based  on 
only  one  design  variable  is  studied.  The  mean  value  of  CH4  concentration  at  the  outlet 
boundary  is  considered  as  the  cost  function.  The  inlet  velocity  is  chosen  as  the  design 
variable.  The  initial  and  boundary  conditions  shown  in  Table  9  are  used  for  the  baseline 
solution.  The  initial  value  for  the  design  variable  is  1 .2  m/s  and  lower  and  upper 
thresholds  (design  constrains)  are  0.3  and  2  m/s  respectively.  The  optimization  is 
performed  by  two  methods:  the  Fletcher-Reeves  conjugate  gradient  method 
(DAKOTA’S  conmin  freg  method)  and  quasi-Newton  method  (DAKOTA’S  optpp  q 
newton  method).  Both  these  methods  are  gradient-based  optimizers  that  are  best  suited 
for  efficient  navigation  to  a  local  minimum  in  the  vicinity  of  the  initial  point. 
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Iteration 


Figure  26.  Convergence  history  of  the  optimization  problem  with  one  design  variable 


—  —  —  Inlet  velocity 

—  •  — Methane  concentration 

wall  temperature 


Figure  27.  Convergence  behavior  of  the  optimization  problem  with  three  design  variables 
using  Fletcher-Reeves  conjugate  gradient  method 
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Figure  26  shows  the  convergence  behavior  of  frcg  and  quasi-Nexton  methods 
for  solving  this  optimization  problem.  As  expected,  inlet  velocity  of  0.3  m/s  is  obtained 
as  the  optimum  value  for  the  design  varible.  The  quasi-Newton  method  shows  the 
quadratic  convergence  rate.  In  the  frcg  method,  7  cost  function  evaluations  and  4 
gradient  evaluations  are  performed.  The  quasi-Newton  used  6  cost  function  evaluations 
and  5  gradient  evaluations. 

The  design  problem  with  three  parameters  including,  inlet  velocity,  inlet 
methane  concentration  and  catalytic  wall  temperature  is  attempted  as  the  second  test 
case.  Since  the  values  of  design  variables  are  spread  across  different  orders  of 
magnitude,  the  normalized  values  are  used  for  the  optimization  cycle.  Constraints 
imposed  in  this  optimization  problem  are  listed  in  Table  12. 


Table  12.  Initial  values  and  design  constraints  used  in  optimization 

Inlet  velocity 

Methane  concentration  at  inlet 

(Normalized  by  multiplying  with  10) 

Catalytic  wall  temperature 

(Normalized  by  dividing  by  1000) 

Lower  bound 

0.3 

0.3 

0.8 

Upper  bound 

1.5 

1.3 

1.2 

Initial  values 

0.7 

0.5 

1.0 

a) 
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Figure  28.  Convergence  behavior  for  the  optimization  problem  with  three  design 
variables  using  the  quasi-Newton  method 
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The  cost  function  is  defined  as  the  mean  value  of  CH4  concentration  at  the  outlet 
boundary.  The  initial  and  boundary  conditions  shown  in  Table  9  are  applied.  Again,  two 
methods  including  Fletcher-Reeves  conjugate  gradient  and  quasi-Newton  are  used  in 
this  problem.  Figures  27  (frcg)  and  28  (quasi-Newton)  show  convergence  behavior  of 
both  these  optimization  algorithms.  In  the  Fletcher-Reeves  conjugate  gradient  method, 
the  cost  function  decreases  from  5.34e-04  to  2.21e-08  for  the  local  optimum  point  of 
(0.3,  0.825,  1.2).  The  cost  function  decreases  from  5.34e-04  to  3.93e-07  using  quasi- 
Newton  method  with  the  local  optimum  point  of  (0.575,  0.825,  1.2). 

As  shown  in  figures  27  and  28,  both  methods  have  a  zig-zag  convergence 
behavior.  The  comparison  of  the  plots  showing  methane  concentartion  along  the 
centerline  of  the  reactor  using  baseline  and  both  optimized  conditions  is  illustrated  in 
figure  29.  The  Fletcher-Reeves  conjugate  gradient  method  is  a  better  optimization 
algorithm  in  this  particular  problem. 


Figure  29.  Comparison  of  methane  concentartion  along  reactor  length 
using  baseline  and  optimized  conditions 


Table  13.  Number  of  solver  and  gradient  calls  for  the  optimization  algorithms 

Method 

Number  of  objective  function 
evaluation  (solver  calls) 

Number  of  objective  gradient 
calculation  (gradient  calls) 

Fletcher-Reeves  conjugate 
gradient  method 

17 

7 

quasi-Newton  algorithm 

15 

13 
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Table  1 3  shows  the  number  of  solver  and  gradient  calls  required  for  both 
optimization  algorithms.  Each  gradient  calculation  step  involves  computing  three 
derivatives  for  this  problem.  Therefore,  the  gradient  calculation  is  the  most  expensive 
part  of  the  optimization  process  and  make  up  about  55%  and  72%  of  the  computational 
cost  of  the  simulation  in  frcg  and  quasi-Newton  algorithms,  respectively.  Since  the 
number  of  design  variables  are  small  in  this  case,  there  is  no  significat  advantage  of 
using  the  adjoint  method  relative  to  the  direct  differentition. 

While  previous  two  caeses  optimized  the  methane  conversion,  the  main  goal  of 
the  reformer  design  is  to  maximize  the  hydrogen  production.  In  some  cases,  although 
methane  is  almost  compeletly  consumed  for  the  given  conditions,  main  products  of  the 
chemistry  are  species  other  than  hydrogen.  For  this  reason,  another  cost  function 
representing  hydrogen  concentration  at  the  outlet  boundary  is  defined  as  following. 

cost  function  =  (59) 


Figure  30.  Convergence  behavior  for  the  optimization  problem  with  three  design 
variables  and  cost  function  representing  hydrogen  concentration  at  the  outlet  boundary 

Based  on  the  previous  experience,  the  Fletcher-Reeves  conjugate  gradient  method 
is  selected  for  this  problem.  Initial  and  boundary  conditions  described  in  Table  9  are  used 
for  the  baseline  case.  Design  variables  are  same  as  the  previous  case  including  inlet 
velocity,  wall  temperature  and  inlet  methane  concentration.  The  initial  values  for  inlet 
velocity  and  wall  temperature  are  same  as  the  previous  test  case.  The  value  of  0.082  is 
chosen  as  the  initial  methane  concentration  at  the  inlet.  Figure  30  shows  the  convergence 
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behavior  for  solving  this  optimization  problem.  Number  of  objective  function 
evaluations  and  gradient  calculations  are  20  and  5,  respectively  to  reach  to  the  local 
minimum  point.  The  comparison  plot  of  hydrogen  concentartion  along  the  centerline  of 
the  reactor  using  baseline  and  optimized  conditions  is  shown  in  figure  3 1 . 


Figure  31.  Comparison  of  hydrogen  concentartion  along  reactor  length  using  baseline  and 

optimized  conditions 


4.4  Iso-Octane  Reforming 

In  this  section,  steam  reforming  results  of  iso-octane  using  model  described  by 
Shi  et  al.  [78]  are  presented.  On-board  hydrogen  production  is  preferred  due  to  lack  of 
infrastructure  for  hydrogen  production,  storage  and  transportation.  Different  methods 
such  as,  partial  oxidation  (PO),  steam  reforming  (SR)  or  autothermal  reforming  (ATR) 
can  be  utilized  for  on-board  production  of  hydrogen  using  commonly  used  liquid  fuels 
such  as  diesel  or  gasoline. 

Partial  oxidation  is  exothermic  reaction  and  does  not  require  external  heat.  Steam 
reforming  is  endothermic  reaction  and  requires  external  heat  to  complete  the  reaction. 
Autothermal  reforming  is  a  combination  of  partial  oxidation  and  autothermal  reforming 
and  is  self-sustaining.  Global  reaction  mechanism  describing  partial  oxidation,  steam 
reforming  and  water-gas  shift  reaction  (WGS)  for  a  general  hydrocarbon  can  be  written 
as  following. 

PO:  CnHm(l)+^02(g)^nC0(g)+^H2(g)  (60.1) 
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SR:  CnHm(l)+nH20(g)^nC0(g)+ 


m  ^ 

— \-n 
2 


H2(g) 


(60.2) 


WGS :  CO(g)  +  H20(g)^C02(g)+H2(g )  (60.3) 

Most  kinetic  models  available  in  literature  employ  methane  (CH4)  as  fuel.  Very 
few  models  have  been  reported  for  heavy  hydrocarbons  such  as,  gasoline  or  diesel  in 
literature  [78,79].  These  models  generally  use  surrogates  to  represent  actual 
hydrocarbons.  For  example, 


(1)  Diesel:  n-hexadecane  (C16H34) 

(2)  Gasoline:  iso-octane  (CsHis) 

Shi  et  al.  [78]  utilized  global  reaction  model  for  n-hexadecane  and  iso-octane  to 
simulate  chemistry  inside  the  reactor.  Thormann  et  al.  [79]  presented  a  detailed  kinetic 
model  for  n-hexadecane.  Model  developed  by  Thormann  et  al.  includes  45  elementary 
reactions,  8  gaseous  species  and  13  adsorbed  species.  Transport  properties  (Lenard-Jones 
parameters)  for  the  hydrocarbons  utilized  in  simulations  can  be  computed  using 
relationships  given  in  a  paper  by  Tee  et  al  [80]. 

A  steam  reforming  model  of  iso-octane  comprising  of  three  global  reactions  and  6 
species  (CsHis,  H2,  CO,  CO2,  ILO,  N2)  is  presented  by  Shi  et  al.  [78].  This  model  is 
implemented  in  the  in-house  multispecies  Navier-Stokes  solver  at  the  SimCenter.  Few 
simulations  have  been  run  using  iso-octane  reforming  model  for  different  values  of  mass 
flow  rates  and  results  are  shown  in  this  report.  Reactions  used  in  iso-octane  reforming 
model  are  given  below.  Details  regarding  reaction  rate  coefficients  and  reaction  rate 
computation  can  be  found  in  the  paper  by  Shi  et  al.  [78]. 


Iso-octane  reforming: 

R\ :  C8//18(/)  +8 H20(g)  ^  8CO(g)+\l H2(g) 

R2:  C0(g)+H20(g)^C02(g)+H2(g) 

R3 :  CgH {S(l)  +  16H 20(g)  <-»  SC02(g)+25H2(g ) 

Figure  32  shows  the  geometry  utilized  in  the  simulation  along  with  different 
boundary  conditions  applied  in  the  problem.  Dimensions  of  the  model  are  given  in  Table 
14.  Operating  conditions  used  in  the  simulation  are  given  in  Table  15. 


Table  14.  Model  dimensions  (Shi  et  al.  [78]) 

Length  (mm) 

Width  (mm) 

Height  (mm) 

200.3 

1.7308 

1.2167 

Table  15.  Operating  Conditions 

Xh20 

Xh20 

XcQ2 

Xh2 

Xc8H18 

Xn2 

T(K) 

P(N/m2) 

55 


303975 


0.005 

0.39 

0.005 

0.005 

0.016 

0.579 

898 

Figure  32.  Reformer  Model  (Shi  et  al.  [78]) 


To  investigate  the  effect  of  fuel  mass  flow  rate  on  temperature  distribution,  few 
cases  have  been  run  with  varying  fuel  velocities  w  hile  keeping  the  rest  of  the  operating 
parameters  constant.  As  seen,  overall  reaction  mechanism  is  endothermic  and  thus, 
reduction  in  temperature  is  evident  in  figures  33  and  34.  However,  extent  of  such 
reduction  is  less  for  the  case  with  higher  velocity  (Figure  33)  by  approximately  30°  K. 
Such  behavior  is  related  to  the  extent  of  chemical  activity  taking  place  inside  the 
reformer.  When  fuel  is  moving  at  slower  velocity  (figure  34),  it  spends  more  time  inside 
the  reformer  and  thus,  has  more  time  to  participate  in  chemical  reactions.  As  overall 
mechanism  is  endothermic,  more  reactions  amount  for  higher  temperature  drop,  which 
justifies  the  trends  shown  in  figure  33  and  34. 


Figure  33.  Temperature  contours  for  V  =  2  m/sec 
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V  =0.25  m/sec 


T(K)  Q  K 


720  74-0  760  780  800  820  840  860  880 


4.5  Fluid-Structure  Interaction  in  SOFC 
4.5.1  Analysis 

The  fluid- structure  interaction  capability  has  been  developed  and  applied  for  the 
SOFC  case  described  in  this  section.  The  geometry  of  the  cell  utilized  in  this  case  is 
shown  in  figure  35.  As  seen,  computational  model  includes  all  relevant  SOFC 
components  including  fuel/air  manifolds,  electrodes/electrolyte  (PEN),  seals  and 
interconnects.  The  number  of  channels  for  both  air  and  fuel  manifolds  is  chosen  to  be 
twelve  to  make  the  geometry  mimic  realistic  planar  SOFC.  Various  geometric 
dimensions  of  the  cell  are  listed  in  Table  16. 


57 


Table  16.  Geometrical  dimensions  of  simplified  SOFC 

Length 

43.0  mm 

Width 

39.0  mm 

Height 

0.23  mm 

Anode  thickness 

0.1  mm 

Cathode  thickness 

0.1  mm 

Seal  thickness 

0.1  mm 

Electrolyte  thickness 

0.1  mm 

Interconnect  thickness 

0.5  mm 

Channel  thickness 

0.5  mm 

No-slip,  adiabatic  wall  boundary  conditions  are  applied  at  the  top  wall,  bottom 
wall  and  side  walls  of  the  computational  geometry  shown  in  figure  35.  Fixed  potential  ( 
(j)  =  0)  boundary  condition  is  applied  at  the  bottom  wall,  while  the  top  wall  is  treated  by 

specifying  average  current  density  ( /  =  iappUed  ).  Inflow  boundary  conditions  with  specified 

mass  flow  rate  and  species  mole  fractions  are  applied  at  both  fuel  and  air  channel  inlets. 
The  temperature  of  the  air  and  fuel  mixture  entering  from  their  respective  inlet  ports  is 
1073  K.  Also,  both  fluids  are  operating  at  atmospheric  pressure.  Specified  back  pressure 
outflow  conditions  are  applied  at  both  air  and  fuel  outlet  ports.  Initial  species  mole 
fractions  and  thermodynamic  conditions  utilized  in  this  simulation  are  given  in  Table  17. 
As  seen  in  Table  17,  partially  reformed  fuel  has  been  utilized  and  thus,  methane 
reforming  and  water  gas  shift  reactions  have  been  considered  inside  the  anode.  Material 
properties  of  different  components  of  SOFC  are  shown  in  Table  18.  Current  density  of 
5500  Am'2  is  applied  at  the  top  wall  of  the  computational  geometry  shown  in  figure  35. 


Table  17.  Mole  fractions  and  thermodynamic  conditions 

^co 

^h2o 

*co, 

** 

XCH, 

Xa, 

** 

T(K) 

P(N  m'2) 

0.029 

0.493 

0.044 

0.263 

0.171 

0.198 

0.802 

1073  K 

101325 

Table  18.  Material  properties  of  various  components  of  SOFC 

Electric  resistivity  of  anode  (flm) 

2.98  xl0‘5exp(- 1392  IT) 

Electric  resistivity  of  cathode  (flw) 

8 . 1 1  x  1  O'5  exp(600  /  T) 

Electric  resistivity  of  interconnect  (flmj 

6.41  xl0‘8 

Ionic  resistivity  of  electrolyte  (flm) 

2.94  x  10"5  exp(10350/  T) 

Thermal  conductivity  of  anode  {w 

40.0 
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/  1 

Thermal  conductivity  of  anode  current  collector  W  m  K 

11.0 

Thermal  conductivity  of  cathode  W 

10.0 

Thermal  conductivity  of  cathode  current  collector  (w 

11.0 

Thermal  conductivity  of  interconnect  (w 

25.0 

Thermal  conductivity  of  electrolyte  [w 

2.0 

Thermal  conductivity  of  seals  [w  rnxK~x 

1.5 

Porosity  of  anode 

0.6 

Porosity  of  cathode 

0.6 

Tortuosity  of  anode 

6.0 

Tortuosity  of  cathode 

6.0 

Porosity  of  anode  current  collector 

0.9 

Porosity  of  cathode  collector 

0.9 

Tortuosity  of  anode  collector 

1.5 

Tortuosity  of  cathode  collector 

1.5 

Two  different  configurations  of  co-flow  and  counter-flow  are  analyzed  in  this 
case.  Figures  36(a)  and  (b)  show  temperature  contours  plotted  over  outer  surface  of  the 
cell  for  co-flow  and  counter-flow  configurations,  respectively.  Figures  also  show  flow 
directions  of  both  air  and  fuel.  As  expected,  in  co-flow  case,  there  is  a  gradual  rise  in 
temperature  as  both  fuel  and  air  move  through  the  flow  domain.  Heat  generated  due  to 
the  electrochemical  reaction  is  the  main  factor  affecting  the  increase  in  temperature.  In 
counter-flow  case,  regions  showing  maximum  temperature  are  present  in  the  middle  of 
the  computational  domain.  Also,  maximum  temperature  found  in  the  co-flow  case  is 
higher  than  the  same  found  in  the  counter-flow  case. 


Temperature  (K):  1080  1100  1120  1140  1180  1180  1200  1 


Figure  36(a).  Co-flow  configuration 


Figure  36(b).  Counter-flow  configuration 


Figure  36.  Temperature  contours  on  outer  surface  on  the  cell 


Temperature  (K):  1080  1100  1120  1140  1160 


Figure  37  shows  polarization  curves  plotted  for  both  co-flow  and  counter-flow 
cases  operating  under  the  same  conditions  described  in  Table  17.  As  expected,  cell 
voltage  reduces  with  increase  in  current  density  due  to  the  effects  of  several 
irreversibilities  present  inside  the  cell.  Both  cases  exhibit  similar  performance  for  low 
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current  densities.  However,  as  current  density  increases,  co-flow  configuration  performs 
better  than  the  counter-flow  configuration. 


Figures  38  (a)  -  (f)  show  mole  fractions  of  different  species  plotted  on  planes 
passing  through  fuel  and  air  manifolds  for  the  co-flow  case.  As  mentioned  earlier,  two 
chemical  reactions  namely,  methane  reforming  and  water-gas  shift  reactions  are 
considered  inside  the  anode  electrode.  Also,  electrochemical  reaction,  which  is 
responsible  for  the  production  of  steam  and  consumption  of  hydrogen  and  oxygen,  affects 
species  distribution  in  the  flowfield.  In  figure  38(a),  there  is  an  overall  reduction  in 
hydrogen  concentration  as  it  moves  through  the  flowfield.  A  region  located  near  bottom 
right  corner  of  the  plane  shows  rise  in  hydrogen  mole  fraction.  This  behavior  is  caused  by 
hydrogen  production  due  to  methane  reforming  reaction.  In  figure  38(b),  gradual  rise  in 
steam  concentration  due  to  electrochemical  reaction  is  evident  as  fuel  moves  from  the 
inlet  to  the  outlet  of  the  manifold.  Methane  mole  fraction  is  plotted  in  figure  38(c).  As 
methane  reforming  is  a  fast  reaction,  most  of  the  methane  can  be  seen  consumed  in  the 
first  half  of  the  flowfield.  Figure  38(d)  shows  oxygen  mole  fraction  plotted  on  a  plane 
extracted  from  the  air  manifold.  As  oxygen  is  a  reactant  of  the  electrochemical  reaction, 
there  is  gradual  reduction  in  its  mole  fraction  as  air  moves  through  the  flowfield. 

Contours  of  carbon  monoxide  (CO)  mole  fraction  plotted  in  figure  38(e)  exhibit  non¬ 
uniformity  in  the  flowfield.  As  CO  acts  as  a  reactant  in  shift  reaction  and  as  a  product  in 
reforming  reaction,  their  combined  effect  produces  non-uniformity  in  the  contours  shown 
in  figure  38(e).  Finally,  contours  of  carbon  dioxide  mole  fraction  are  plotted  in  figure 
38(f).  As  the  only  reaction  involving  carbon  dioxide  (as  a  product)  is  a  shift  reaction, 
gradual  rise  in  carbon  dioxide  concentration  is  evident  in  figure  38(f). 
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Figures  39(a)  -  (f)  show  stress  contours  plotted  on  planes  extracted  through 
different  solid  and  porous  components  of  the  cell  for  co-flow  configuration.  Figures  39(a) 
and  (b)  show  contours  of  maximum  tresca  equivalent  stress  (MTES)  plotted  on 
streamwise  planes  passing  through  seals  located  on  both  fuel  and  air  sides,  respectively. 
As  seen,  maximum  stress  is  present  near  fuel  inlet  port  in  both  figures.  Figure  39(c)  -  (e) 
show  contours  of  mean  principal  stress  (MPS)  plotted  on  planes  passing  through  the 
anode ,  electrolyte  and  cathode,  respectively.  Some  characteristics  of  these  contours  such 
as  location  of  the  maximum  stress  region  are  similar  in  all  three  plots.  This  region  is 
located  near  fuel  inlet  port  in  the  computational  domain.  Figure  39(f)  shows  contours  of 
MTES  plotted  on  a  vertical  plane  passing  through  the  interconnect  and  inlet  ports.  As 
expected,  maximum  stress  is  found  near  fuel  inlet  port.  Regions  near  air  inlet  port  also 
ind.cale  high  values  of  MTES.  Overall,  maximum  stress  values  are  found  in  a  plane 
extracted  through  the  anode  electrode  in  figure  39(c). 
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(a).  Anode  seal 


(b).  Cathode  seal 


(c).  Anode 


(d).  Electrolyte 


(e).  Cathode 


MTES(MPa):  30  50  70  90  110  130 


a 

EJ 


(Q.  Interconnect 


Figure  39.  Stress  distribution  inside  different  components  for  co-flow  configuration 


Figures  40(a)  -  (f)  show  stress  contours  plotted  on  planes  extracted  through 
different  components  of  the  cell  for  counter-flow  configuration.  Even  though  air  is 
flowing  in  the  opposite  direction  in  this  case  compared  to  the  previous  case  (figure  39), 
stress  contours  in  different  components  show  similar  characteristics  in  both  cases. 
Regions  with  maximum  stress  are  located  near  fuel  inlet  port  in  different  cell 
components.  Overall,  stress  values  for  the  counter-flow  case  are  smaller  than  the  co-flow 
case. 


- 


(b).  Cathode  seal 
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(c).  Anode 


(d).  Electrolyte 


(e).  Cathode 


□ 


ipT- 


(f).  Interconnect 


Figure  40.  Stress  distribution  inside  different  components  for  counter-flow  configuration 


4.5.2  Stress  Sensitivity 

As  mentioned  earlier,  both  fuel  cell  and  structures  code  are  capable  of  performing 
sensitivity  analysis.  To  demonstrate  this  capability,  stress  sensitivity  contours  are  plotted 
on  planes  extracted  through  different  components  of  the  cell  in  figures  41(a)  -  (e).  The 
co-flow  configuration  is  utilized  in  figure  41 .  Design  variable  in  this  study  is  cathode 
porosity.  The  method  utilized  to  compute  sensitivity  derivatives  is  direct  differentiation  in 
both  fuel  cell  and  structures  code.  Structures  code  requires  values  of  flowfield  variables 
and  sensitivities  of  flowfield  variables  from  the  fuel  cell  code  to  compute  stress 
sensitivities.  The  characteristics  of  stress  sensitivity  contours  in  figure  41  are  similar  to 
those  shown  for  the  stress  contours  in  figure  39;  especially  regions  with  highest 
sensitivity  values  are  located  near  fuel  inlet  port  for  all  components. 
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(c).  Anode 


(d).  Electrolyte 


(e).  Cathode 


Figure  41 .  Stress  sensitivities  with  respect  to  the  cathode  porosity  inside  different  components  for  co-flow 
_ configuration _ 


4.6  CAD  Integration  and  Mesh  Movement 

As  described  in  section  3.4.2,  a  framework  has  been  developed  for  integrating 
CAPRI  [67]  with  the  SimCenter  geometry  libraries  so  that  CAD-based  design  variables 
can  be  utilized.  In  this  section,  three  different  examples  demonstrating  use  of  CAPRI 
interface  and  modified  linear  elasticity  mesh  movement  solver  are  provided. 

First  example  involves  modifying  the  shape  of  the  air  channel  for  the  tubular 
SOFC  depicted  in  figure  42.  Here,  the  geometry  has  been  constructed  in  a  CAD  program 
(Solid  Works)  and  the  air  channel  is  described  as  a  simple  ellipse,  where  the  dimensions 
of  the  vertical  and  horizontal  axes  are  exposed  as  design  variables.  The  original  design,  as 
shown  in  the  upper  left  portion  of  figure  43.  has  an  air  channel  w  ith  a  circular  cross 
section.  Sensitivity  derivatives  describing  the  changes  in  the  surface  coordinates  with 
each  design  variable  are  shown  in  the  right  portion  of  this  figure,  whereas  the  result  of  a 
simple  design  problem  is  shown  in  the  lower  left.  As  seen  in  the  final  design,  the  circular 
cross  section  has  become  elliptical  after  one  design  cycle  and  the  fidelity  of  the  mesh  is 
maintained  after  repositioning  the  coordinates  of  the  air  channel. 
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Air  Channel 


Figure  42.  Shape  design  of  tubular  SOFC 


Figure  43.  Shape  design  and  sensitivity  contours  of  tubular  SOFC 


Second  example  [81]  is  demonstrated  in  figure  44  showing  capability  to  produce 
quadratic  grids  using  CAPRI  and  linear  elasticity  solver.  Shape  shown  in  figure  44  is  an 
analytically  defined  three-dimensional  body  of  revolution  [82],  Figure  44(a)  shows  the 
high-order  mesh  generated  without  getting  original  definition  of  the  CAD  model.  The 
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linear  surface  representation  shown  in  figure  44(a)  is  clearly  inaccurate  and  using  it  for 
higher-order  schemes  would  not  allow  accurate  definition  for  the  geometry.  In  contrast, 
figure  44(b)  shows  a  surface  mesh  for  the  same  geometry  where  quadratic  elements  are 
used.  After  adding  additional  surface  quadrature  points  using  the  CAPRI  interface  to  the 
elements,  the  fidelity  of  the  surface  is  clearly  improved. 


(a)  Linear  surface  (b)  Quadratic  surface 


Figure  44.  Surface  representation  for  an  analytic  3D  body  of  revolution  [81,82] 

Third  example  illustrates  the  mesh  movement  strategy  for  a  NACA  4412  mesh 
shown  in  figure  45(a)  [81].  The  mesh  is  generated  with  a  viscous  spacing  of  le-03  normal 
to  the  wall.  Figure  45(b)  plots  contours  of  the  perturbation  magnitude  obtained  using 
linear  elasticity  solver  described  in  Section  3.5.  As  expected,  most  perturbations  are 
concentrated  on  the  airfoil  upper  surface  because  of  the  presence  of  relatively  larger 
curvatures,  and  the  perturbation  magnitude  quickly  decays  as  the  distance  to  the  airfoil 
increases.  Figures  45(c)  and  (d)  depict  a  close-up  view  of  the  viscous  layer  near  0.4  chord 
of  length  on  the  upper  airfoil  surface  for  linear  and  higher-order  mesh,  respectively.  As 
shown  in  figure  45(d),  the  mesh  movement  strategy  effectively  produces  sufficient 
deformations  for  the  interior  mesh  points  and  quadrature  points  to  prevent  negative 
volumes  or  Jacobians.  The  mesh  movement  strategy  employed  in  the  current  w'ork  is 
capable  of  obtaining  valid  finite-element  meshes  with  high-aspect  ratio  elements  for  high 
Reynolds  number  flows. 


(a)  NACA  4412  mesh  (b)  Displacements  obtained  using  LE  solver 
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if 


(c)  Linear  mesh 


(d)  Higher-order  mesh 


Figure  45.  Mesh  movement  strategy  for  higher-order  grids  using  CAPRI  and  linear 

elasticity  solver  [8 1  ] 


5.  Conclusions 

An  in-house  solver  capable  of  solving  multispecies  Navier-Stokes  equations  has 
been  developed.  The  solver  is  implicit,  unstructured  and  implemented  in  parallel  for 
computational  efficiency.  Two  applications  explored  in  this  project  include  planar  SOFC 
and  catalytic  reactor.  Results  obtained  using  numerical  models  of  SOFC  as  well  as 
catalytic  combustion  of  methane  are  validated  by  performing  comparisons  with  the 
experimental  results  from  the  literature.  Effects  of  inlet  mass  flow  rate  conditions  on 
performance  of  both  SOFC  and  reformer  are  investigated.  The  surface  chemistry, 
heterogeneous  combustion  and  coverages  are  computed  by  linking  the  solver  with 
Cantera.  Parametric  studies  are  carried  out  for  several  design  parameters  affecting  reactor 
efficiency  for  fuel  reforming.  The  code  is  linked  to  DAKOTA  for  computational  design 
and  optimization  purposes.  Based  on  the  presented  results,  the  following  conclusions  may 
be  drawn. 

The  parameter  study  shows  that 

o  Rhodium  shows  better  performance  for  partial  oxidation  of  methane 
(conversion  of  90%)  than  platinum  (conversion  of  77%). 

o  Conversion  of  methane  is  predicted  to  decrease  with  increasing  inlet 
velocity  and  Reynolds  number. 

o  Conversion  of  methane  and  hydrogen  production  increase  with 
increase  in  catalytic  wall  temperature. 

Sensitivity  analysis  for  the  methane  reformer  shows  that  the  methane 
concentration  is  an  important  parameter  affecting  reformer  performance. 
Sensitivity  derivatives  are  used  for  computational  design  based  on  gradient- 
based  optimization  algorithm.  Computation  of  sensitivity  derivatives  using 
direct  differentiation  and  adjoint  method  reduced  the  run  time  of  the 
optimization  process  up  to  38%  for  the  Fletcher-Reeves  conjugate  gradient 
method  and  54%  for  the  quasi-Newton  algorithm. 

A  capability  to  perform  thermo-mechanical  analysis  of  different  components  of 
SOFCs  has  been  developed  by  coupling  the  multispecies  solver  with  the  structures  code. 
Results  obtained  using  this  capability  indicate  that  the  main  factors  affecting  the  stress 
distribution  are  temperature  gradients  and  mismatch  of  coefficients  of  thermal  expansion 
between  different  components  of  the  cell.  Also,  capability  to  perform  thermo-mechanical 
sensitivity  analysis  has  been  developed  using  direct  differentiation  method.  Using  this 
capability,  sensitivity  derivatives  of  stress  has  been  computed  for  planar  SOFC. 
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Capability  to  perform  shape  design  and  create  curved  higher  order  elements  using 
original  definition  of  CAD  model  is  developed  using  CAPRI  (developed  at  MIT)  and 
high-fidelity  mesh  deformation  algorithms.  This  method  is  successfully  applied  to  create 
higher  order  curved  elements  for  viscous  grid  of  N  AC  A  4412. 
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Nomenclature 


Symbols 

Name 

Unit 

B 

Permeability 

m2 

E, 

Total  energy 

Jm'3 

f 

Cost  function 

cost  function  depended 

H 

Enthalpy 

J  kg'1 

i 

Current  density 

Am'2 

K 

Exchange  current  density 

Am'2 

j 

Mass  flux  vector 

kg  m'2  s'1 

• 

Mass  flow  rate 

kgs'1 

m 

kg  kmol'1 

M 

Molecular  weight 

ns 

Number  of  species 

- 

P 

Pressure 

Nm'2 

Q 

Solution  vector 

solution  variable  dependent 

q 

Heat  flux 

J  m'2  s'1 

T 

Temperature 

K 

u 

x-velocity  component 

m  s'1 

V 

y-velocity  component 

m  s'1 

w 

z- velocity  component 

m  s'1 

x, 

Mole  fraction  of  i th  species 

- 

Greek  Symbols 

P 

Design  variable  vector 

design  variable  dependent 
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X 

Grid  vector 

m 

p 

Mass  concentration 

kg  m'3 

p 

Molecular  viscosity 

kg  m  s'1 

<p 

Electric  potential 

volt 

9 

Activation  polarization 

volt 

e 

Porosity 

- 

K 

Totruosity 

- 

V 

Control  volume 

m3 

r 

Viscous  flux 

kg  m'1  s'2 

Constants 

F 

Faraday  constant 

96484.56  A  s  mol'1 

K 

Universal  gas  Constant 

8314.4  Jkmol'1  K'1 

Indices 

a 

Anode 

c 

Cathode 

eff 

Effective 

i,j 

Chemical  species 
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1.  Executive  summary 

In  order  to  supply  hydrogen  or  syngas  for  fuel  cell-based  auxiliary  power  units,  onboard  fuel 
processing  technology  has  been  paid  significant  attention  in  energy  and  fuel  cell  society.  For 
providing  hydrogen  rich  fuels  to  fuel  cells,  converting  hydrocarbon  fuels  into  syngas  through 
reforming/shift  reactions  has  relatively  less  technical  obstacles  regarding  the  storage  of 
hydrocarbon  fuels  compared  to  hydrogen  gas  fuel. 

This  project  focused  on  onboard  fuel  processing  of  commercial  Jet-A  (equivalent  to  JP-8 
regarding  sulfur  and  hydrocarbon  content)  fuel  to  produce  hydrogen  and  syngas  for  fuel  cell 
auxiliary  power  units.  Jet-A  fuel  was  studied  because  it  is  the  logistic  fuel  (equivalent  to  JP-8) 
commonly  used  for  civilian  airplanes  and  military  heavy  duty  trucks.  The  research  team 
developed  a  new  catalyst,  using  low  cost  materials,  for  desulfurization  from  Jet-A  fuel  at  room 
temperature  and  ambient  pressure,  which  is  innovative  and  cost-effective.  Ultra-deep  adsorptive 
desulfurization  has  been  achieved  for  Jet-A  fuel  from  over  1,000  ppmw  to  below  50  ppmw.  The 
second  part  of  this  work  focused  on  autothermal  reforming  of  desulfurized  jet-A  fuel. 
Autothermal  reforming  catalyst  and  monolith  reaction  section  have  been  developed.  The 
experimental  tests  for  syngas  production  were  conducted  firstly  using  rc-dodecane  (as  a  surrogate 
of  Jet-A  fuel),  and  then  using  desulfurized  Jet-A  fuel.  A  reactor  accommodating  both  reformer 
and  water-shift  reactions  for  desulfurized  Jet-A  fuel  has  been  designed  and  operated. 

For  the  adsorptive  desulfurization  of  Jet-A  fuel,  a  novel  cost-effective  NiO-CeCVALCL-SiCL 
adsorbent  was  proposed  and  prepared  in-house  for  experimental  tests.  Details  of  the 
compositions  of  catalyst  and  procedures  of  making  and  processing  the  catalyst  are  provided  in 
this  report.  The  sulfur  adsorption  kinetic  characteristic  and  isotherm  at  equilibrium  were  studied 
in  batch  tests.  The  maximum  desulfurization  efficiency  could  reach  up  to  98%.  For  Jet-A  fuel  in 
a  total  sulfur  concentration  of  1037  ppmw,  the  lowest  achieved  sulfur  concentration  at  room 
temperature  in  the  treated  fuel  was  22.13  ppmw.  This  is  a  significant  achievement  regarding  the 
desulfurization  efficiency,  especially  at  room  temperature  and  using  commercial  Jet-A  fuel  with 
no  pretreatment.  The  dynamic  desulfurization  performance  of  the  adsorbent  was  investigated  in 
fixed-bed  reactor  tests.  The  flow  rate  of  fuel,  size  of  adsorbent  particles,  and  dimensions  of  an 
adsorbent-packed  fixed-bed  were  optimized  to  obtain  a  high  sulfur  adsorption  capacity  per  unit 
mass  of  adsorbent.  At  a  breakthrough  sulfur  concentration  of  10  ppmw  a  very  high  sulfur 
adsorption  capacity  of  0.633  mg  S/g  adsorbent  was  achieved.  At  a  mean  sulfur  concentration  of 
30  ppmw  in  the  treated  accumulated  fuel,  the  best  capacity  achieved  is  1.98  mg  S/g  adsorbent. 
The  scaling-up  strategies  for  the  fixed-bed  reactor  and  the  method  for  adsorbent  regeneration  are 
also  investigated.  Finally,  preliminary  tests  of  adsorbent  regeneration  showed  that  the  first-time 
regenerated  adsorbent  could  recover  83%  of  the  sulfur  removal  capacity  compared  to  a  fresh 
adsorbent,  and  a  second  time  regeneration  could  recover  50.4%  of  that  of  fresh  adsorbent. 

For  the  reforming  of  Jet-A  fuel,  autothermal  reforming  (ATR)  method  was  employed  and  a 
bimetallic  NiO-Rh  catalyst  with  promoters  of  Ce,  K,  and  La  were  synthesized  for  the  ATR 
reactions.  A  lab-scale  2.5  kWt  autothermal  reforming  system  including  the  reformer  and  balance- 
of-plant  was  designed,  fabricated,  integrated  and  tested.  The  reforming  system  performance  at 
various  operation  conditions  was  compared.  Reformer  operation  temperature,  steam  to  carbon 
ratio  and  oxygen  to  carbon  ratio,  as  well  as  pre-heating  temperatures  for  fuel,  air,  and  steam  were 
optimized  for  better  system  energy  conversion  efficiency,  H2  selectivity,  and  COx  selectivity. 
For  n-dodecane  the  energy  conversion  efficiency  was  83.5%  while  for  desulfurized  commercial 
Jet-A  fuel  the  efficiency  reduced  to  77.0%. 
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2.  Introduction 

High  energy  efficiency  and  energy  density,  together  with  rapid  refuel  capability,  render  fuel 
cells  being  highly  attractive  for  portable  power  generation,  serving  as  auxiliary/backup  power 
units  (APUs).  Fuel  cell  APUs  need  hydrogen  or  syngas  as  the  energy  source  to  produce  electrical 
power,  however,  there  is  no  existing  infrastructure  for  hydrogen  production  and  storage.  To 
provide  hydrogen  for  fuel  cells  in  military  and  civilian  APUs,  jet  fuels  are  widely  considered  as 
an  excellent  choice  since  they  are  logistic  fuels  for  both  civilian  and  military  transportation 
vehicles.  Integrated  micro  fuel  processors  in  combination  with  solid  oxide  fuel  cell  (SOFC) 
stacks  using  jet  fuels  have  been  viewed  as  achievable  technologies  for  power  generation. 

There  are  various  organic  sulfur  compounds  contained  in  jet  fuels.  As  regulated  by  the  US 
Environmental  Protection  Agency  (EPA)  in  2010,  the  highest  sulfur  concentration  allowed  in  jet 
fuels  is  3000  ppmw.  The  major  sulfur  compounds  in  jet  fuels  are  2,3-DMBT  (2,3- 
dimethylbenzothiophene),  2,3,7-TMBT  (2,3,7-trimethylbenzothiophene),  2,3,5-TMBT  (2,3,5- 
trimethylbenzothiophene),  and  2,3,6-TMBT  (2,3,6-trimethylbenzothiophene).  The  sulfur 
contents  must  be  removed  before  jet  fuels  can  be  reformed  to  produce  hydrogen  or  syngas  for 
fuel  cells.  This  is  because  the  presence  of  high  level  sulfur  compounds  is  harmful  to  the  catalyst 
for  the  reforming  reaction,  and  also  because  sulfur  is  poisonous  to  the  electrodes  and  catalysts  in 
fuel  cells.  Therefore,  for  jet  fuels  being  supplied  to  reforming  and  water-shift  for  producing 
hydrogen  rich  fuels,  there  are  three  critical  research  subjects  to  be  addressed:  1)  catalytic 
desulfurization  of  the  fuel;  2)  autothermal  catalytic  reforming  with  precisely  controlled  oxygen, 
steam,  and  fuel  at  a  temperature;  3)  system  integration  and  thermal  management  for  reforming 
and  water-shift  to  be  included  in  one  reactor. 

Conventional  desulfurization  technology  employed  in  refineries  is  the  hydrodesulfurization 
(HDS)  method.  This  technology,  working  at  300-400  °C  and  40-50  bars,  requires  heavy 
expenditure  in  both  capital  and  operation.  It  also  needs  large  volume  reactors,  making  it 
inconvenient  for  portable  fuel  cell  applications.  Another  drawback  of  HDS  is  that  the  octane 
number  is  significantly  reduced,  due  to  the  saturation  of  alkenes  and  arenes  by  hydrogen  at  high 
temperatures  and  pressures.  Desulfurization  of  commercial  jet  fuel  at  room  temperature  and 
atmospheric  pressure  is  desirable  for  practicality  and  simplicity  of  the  system.  Selective 
adsorptive  desulfurization  or  (selective  adsorption  for  removing  sulfur — SARS)  of  real  Jet-A 
fuel  with  an  initial  total  sulfur  concentration  of  over  1000  ppmw  at  room  temperature  and 
atmospheric  pressure  is  the  objective  of  our  present  study.  Using  commercial  Jet-A  fuel  instead 
of  model  fuel  is  a  big  challenge  in  the  adsorptive  desulfurization  study,  because  aromatics  and 
olefins  in  the  real  fuel  have  a  strong  inhibiting  effect  on  the  adsorptive  desulfurization 
performance,  especially  at  low  temperatures  [1]. 

Desulfurization  methods  suitable  for  portable  applications  without  reducing  the  fuel  octane 
number  have  been  studied  in  recent  years,  and  two  of  the  promising  processes  were 
pervaporation  [2-4]  and  selective  adsorption  [5-8].  The  latter  is  considered  the  most  promising 
method  [1].  Various  adsorbents,  based  on  transition  metals  and  supported  materials  including 
metals,  carbon,  and  zeolites,  have  been  proposed  for  sulfur  adsorption  performances  [9-13]. 
Nickel  and  copper  supported  on  activated  carbon  and  zeolite  exhibited  effective  desulfurization 
effects  for  jet  fuels  [14,15].  Muzic  et  al.  [16]  examined  the  adsorptive  desulfurization  of  three 
commercial  activated  carbon  adsorbents  on  diesel  fuel.  Dastanian  and  his  coworkers  [17]  studied 
the  desulfurization  of  gasoline  containing  140  ppmw  total  sulfur  by  using  a  nanoporous  Ni 
loaded,  Y-type  zeolite  at  ambient  conditions.  The  adsorption  capacity  of  Ni-Y  zeolite  was  0.84 
and  2.31  mg  S/g  adsorbent,  corresponding  to  the  fuel  residual  sulfur  concentrations  of  10  and  30 
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ppmw,  respectively,  in  a  batch  test.  Montazerolghaem  et  al.  [18]  investigated  the  adsorptive 
sulfur  removal  from  a  modeled  gasoline  by  Na-Y  zeolite  and  Ce-Y  zeolite  at  room  temperature. 
Their  sulfur  model  solution  consisted  of  iso-octane  and  thiophene  that  contained  about  116 
ppmw  total  sulfur  compounds  initially,  and  after  the  desulfurization,  the  total  sulfur 
concentration  was  about  20  ppmw.  The  Ce-Y  zeolite  was  reported  to  have  a  better 
desulfurization  performance  than  the  Na-Y  zeolite,  and  the  capacity  of  Ce-Y  zeolite  was  around 
0.68  mg  S/g  adsorbent  at  the  adsorbent  to  modeled  gasoline  ratio  of  1  g/10  ml  fuel.  Nevertheless, 
these  mentioned  catalysts  and  supporter  materials  (such  as  activated  carbon  and  zeolite)  are 
generally  expensive,  and  search  for  low  cost  catalyst  for  practical  application  is  focused  in  this 
project. 

The  ideal  adsorbent  for  the  SARS  process  at  room  temperature  and  atmospheric  pressure  must 
simultaneously  satisfy  merits  of  strong  Lewis  surface  acidity,  big  specific  surface  area,  strong 
ionic  polarity,  perfect  redox  properties,  and  chemical  and  thermal  stabilities.  In  our  works  [12], 
using  the  advantages  of  the  SARS  mechanisms  of  p-complexation,  direct  S-M  interaction,  active 
sites  for  adsorption,  and  the  promotional  effect  of  adsorbent  material  structure  corresponding  to 
the  novel  geometric  effect  and  the  meso-porous  structure,  we  developed  a  promising  adsorbent, 
Ni-Ce/ALCb-SiCL,  for  the  ultra-deep  desulfurization  of  jet  fuels  with  a  high  sulfur  content  at 
room  temperature.  The  adsorbent  to  be  investigated  in  this  work  has  been  proven  to  have  good 
sulfur  compounds  selectivity,  reasonable  adsorptive  capacity,  and  the  ability  to  decrease  the  total 
sulfur  concentration  to  below  50  ppmw.  Very  importantly,  the  cost  of  the  catalyst  materials  is 
low  which  important  for  large  scale  application. 

The  three  main  approaches  for  reforming  of  hydrocarbon  fuels  are  steam  reforming  (SR), 
partial  oxidation  (POX),  and  autothermal  reforming  (ATR).  There  has  been  a  number  of  research 
works  studied  the  reforming  of  different  fuels,  such  as  ethanol,  gasoline,  diesel,  and  jet  fuels 
employing  one  of  the  above-mentioned  approaches.  It  is  understood  from  the  literature  that  SR 
has  high  H2  concentration  in  reformate,  but  it  requires  lots  of  external  heat  and  the  system  is 
usually  bulky  and  heavy.  POX  has  a  compact  system  and  the  reaction  is  exothermic,  however, 
the  H2/CO  ratio  in  reformate  is  relatively  low.  ATR  is  a  combination  of  POX  followed  by  SR 
and  it  has  the  most  suitable  characteristics  for  onboard  fuel  reforming.  An  ATR  system  has 
favorable  H2/CO  ratio  in  reformate,  less  coke  formation  tendency,  no  requirement  of  external 
heat  source,  relatively  compact  system  size  and  weight,  and  rapid  startup  and  dynamic  responses. 
Because  of  the  advantages  of  ATR,  it  was  selected  as  the  approach  of  reforming  in  the  present 
study.  Both  reforming  catalyst  and  operation  conditions  can  significantly  influence  the  ATR 
reaction  in  respect  to  H2  and  CO  concentrations,  hydrogen  and  COx  selectivity,  as  well  as  system 
energy  conversion  efficiency  [19],  Ni-based  catalysts  for  fuel  reforming  have  been  widely  used 
because  of  the  high  activity  and  low  cost.  However,  Ni  catalysts  have  inherent  challenges  such  as 
sulfur  poisoning  and  coke  formation  [20].  Noble  metals  were  proven  to  be  effective  in  fuel 
reforming  reactions,  but  the  high  price  of  noble  metals  is  a  disadvantage  [21,22].  Thus  bimetallic 
catalysts  (including  nickel  and  noble  metal)  have  been  proposed  and  investigated  in  some  works 
for  fuel  reforming.  Effects  of  ATR  operating  conditions  such  as  steam  to  carbon  ratio  (S/C)  and 
oxygen  to  carbon  ratio  (O2/C)  were  investigated  both  numerically  and  experimentally  for 
different  reformers,  and  it  is  believed  that  the  reasonable  working  ranges  of  S/C  and  O2/C  are 
1.25-2.5  and  0.35-0.5,  respectively.  A  lab-scale  2.5  kWt  autothermal  reforming  system  with  a 
new  reformer  design  and  novel  catalyst  was  experimentally  studied.  The  used  bimetallic  catalyst 
is  NiO-Rh  which  is  mixed  with  promoters  of  CeC>2,  K2O,  and  La2C>3,  and  was  prepared  in-house 
for  experimental  test.  Using  n-dodecane  as  a  surrogate  of  Jet-A  fuel,  the  effects  of  operating 
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conditions  such  as  reformer  temperature,  S/C,  and  O2/C  were  experimentally  investigated  in 
respect  to  hydrogen  selectivity,  COx  selectivity,  and  the  energy  conversion  efficiency.  Coke 
formation  was  suppressed  by  reducing  the  pre-heating  temperature  of  the  fuel.  After  the 
optimized  operating  conditions  were  determined,  desulfurized  commercial  Jet-A  fuel  was  tested 
in  the  ATR  system, 

Because  of  the  complexity  of  carbon  content  in  jet- A  fuel,  surrogate  of  Jet-A  fuel,  n- 
dodecane,  was  firstly  used  in  laboratory  experiments  to  investigate  the  reforming  characteristics. 
Real  desulfurized  commercial  Jet-A  fuel  was  thereafter  tested  at  the  optimized  reforming 
operation  conditions. 

3.  A  new  selective  adsorption  catalyst  for  Jet  fuel  desulfurization  at  room  temperatures 

The  main  objectives  of  the  present  work  include:  1)  addressing  the  challenges  of  removing 
sulfur  compounds  in  commercial  jet  fuels  to  an  acceptable  level  using  low  cost  catalyst;  2) 
reforming  the  low-sulfur  jet  fuels  into  hydrogen  rich  syngas  for  use  in  SOFC  APUs.  Jet-A  and 
JP-8  are  the  two  most  commonly  used  jet  fuels  in  civilian  and  military  applications.  Since  JP-8 
fuel  is  for  military  use  and  restricted  for  public  to  purchase,  Jet-A  fuel  was  selected  as  the  test 
fuel  in  this  study.  The  hydrocarbon  and  sulfur  content  of  jet- A  and  JP-8  are  believed  the  same. 

3.1  Jet  fuel  analysis 

Jet  fuels  are  commonly  used  in  heavy-duty  trucks,  aircraft,  and  ships,  and  they  are  also  the 
preferred  logistic  fuels  in  military  applications  because  of  their  higher  efficiency  and  power 
density  as  well  as  lower  flammability  than  gasoline.  JP-8  was  selected  as  the  exclusive  battlefield 
fuel  by  the  Department  of  Defense  and  North  Atlantic  Treaty  Organization  [23]  due  to  its  high 
flash  point  of  38  °C  and  good  low-temperature  operations  [24].  The  specifications  for  JP-8  are 
similar  to  Jet-A  except  that  JP-8  has  required  additives  for  anti-icing,  a  corrosion  inhibitor, 
lubricity  improver,  and  anti-static,  as  well  as  antioxidant  and  metal  deactivators  [25].  Jet-A  is  a 
kerosene  type  of  fuel  and  is  widely  used  as  fuel  for  civilian  airliners  [26].  A  better  understanding 
of  diesel  and  jet  fuel  properties  will  be  helpful  to  the  development  of  an  onboard  fuel  processing 
system  for  the  fuel  cell  APUs. 

Jet  fuels  consist  of  thousands  of  chemicals,  mainly  hydrocarbons,  as  well  as  functional 
additives.  Kerosene  types  Jet-A  and  JP-8  fuels  have  approximate  carbon  number  distributions 
between  8  and  16,  and  H/C  molar  ratios  from  1.6  to  2.0.  The  average  molecular  formula  of  JP-8 
fuel  is  Ci,H2i  [27]  and  the  average  molecular  formula  of  Jet-A  fuel  is  C11.6H22.3  [28],  Aromatics 
in  jet  fuels  are  limited  to  no  more  than  25%  by  volume  [29],  As  listed  in  reference  [20].  typical 
jet  fuels  in  the  U.S.  has  an  average  71%  volume  of  paraffins,  19%  volume  of  aromatics,  6.2% 
volume  of  naphthalenes  and  3.5%  volume  of  olefins. 

Various  sulfur  compounds  are  contained  in  commercial  jet  fuels.  The  major  sulfur  compounds 
in  jet  fuels  are  2,3-DMBT  (2,3-dimethylbenzothiophene),  2,3,7-TMBT  (2,3,7- 

trimethylbenzo thiophene),  2,3,5-TMBT  (2,3,5-trimethylbenzothiophene),  and  2,3,6-TMBT 
(2,3,6-trimethylbenzothiophene).  Desulfurization  will  remove  these  compounds  as  much  as 
possible. 

Sulfur  compounds  contained  in  both  diesel  and  jet  fuels  are  unwanted  because  they  are 
poisonous  to  catalysts  and  electrodes  in  fuel  cells.  In  the  U.S.,  the  regulations  for  sulfur  content 
in  diesel  and  jet  fuels  set  by  the  Environmental  Protection  Agency  (EPA)  have  become  more  and 
more  strict.  Table  1  shows  the  regulations  for  sulfur  content  from  1993  to  2010  [30].  Although 
the  sulfur  concentration  regulation  set  by  EPA  in  2010  is  less  than  3,000  ppmw,  the  average 
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value  of  sulfur  concentration  for  jet  fuel  reported  in  the  literature  is  714  ppmw.  However,  there  is 
a  standard  deviation  of  414  ppmw,  which  reflects  that  the  variations  of  sulfur  concentrations  in 
different  fuel  sources  are  quite  large  [31-36].  The  jet-A  fuel  in  this  work  was  purchased  from 
Million  Air,  Tucson — a  local  fuel  supplier  at  the  airport.  The  total  sulfur  concentration  in  the  fuel 
was  measured  by  a  Thermo  TS  3000  total  sulfur  analyzer,  which  has  a  working  range  of  0.02- 
5000  ppmw  and  an  uncertainty  of  less  than  5%  of  the  measured  value.  Our  measurement  to  the 
commercial  jet-A  fuel  showed  a  content  of  1 140  ppmw  of  sulfur. 


Table  7  U.S.  EPA  ’s  sulfw -  regulations  for  diesel  andjet £uels  [ 30] 


Category 

Year 

1993 

2006 

2010 

Highway  diesel  (ppmw) 

500 

15 

15 

Non-road  diesel  (ppmw) 

5000 

500 

15 

Jet  fuel  (ppmw) 

3000 

3000  max 

<3000 

Surrogates  of  diesel  and  jet  fuels  are  usually  used  in  laboratory  research  and  experiments  to 
better  calculate  and  control  the  reforming  processes  for  optimization.  Surrogates  are  typically 
divided  into  two  types.  Physical  surrogates  are  mixtures  that  have  similar  physical  properties 
(e.g.,  heating  value,  density,  viscosity,  specific  heat  capacity,  flash  point,  etc.)  to  real  fuels. 
Chemical  surrogates  are  mixtures  that  have  generally  the  same  chemical-class  compositions  (e.g., 
aromatics,  naphthenes,  olefins,  alkanes,  etc.)  and  average  molecular  weight  as  real  fuels  [37,38]. 
In  the  fuel  reforming  process,  surrogates  are  expected  to  reproduce  most  of  the  fuel  chemical 
characteristics,  as  well  as  some  important  physical  properties  in  heat  and  mass  transfer  [39]  The 
lower  heating  value  (LHV)  is  one  most  important  property  of  fuels  since  it  determines  the 
reforming  efficiency.  The  quantity  of  LHV  is  determined  by  subtracting  the  heat  of  water 
vaporization  from  the  higher  heating  value  (HHV),  which  is  the  thermodynamic  heat  of 
combustion.  A  number  of  recent  studies  on  surrogates  of  diesel  and  kerosene  type  jet  fuels  for 
autothermal  reforming  have  been  reported  [40-45]  n-dodecane  (>99%)  was  proved  to  be  a  good 
surrogate  of  Jet-A  [39]  to  study  the  autothermal  reforming  in  terms  of  fuel  conversion,  hydrogen 
yield,  reactor  temperature  profile  and  reforming  efficiency.  As  a  surrogate  fuel,  n-dodecane  has 
the  chemical  formula  of  C12H26  with  hydrogen  content  15.28  wt.%.  And  the  molecular  weight 
was  170.3  g/mol,  and  LHV  was  44.14  MJ/kg  (LHV  of  Jet-A  fuel  is  43.26  MJ/kg). 

3.2  Proposed  catalyst  and  the  procedures  of  preparation 

3.2.1  Design  of  catalyst 

Two  most  popular  types  of  adsorptive  desulfurization  mechanisms  are  7r-complexation  [46-48] 
and  direct  Sulfur-Metal  (S-M)  interaction.  The  71-complexation  adsorbents  particularly  the  Y- 
zeolites  exhibit  high  sulfur-adsorption  capacity,  but  show  low  selectivity  for  sulfur  compounds 
as  the  result  of  competitive  adsorption  of  aromatic  compounds.  Meanwhile,  according  to  the 
adsorption  mechanism  of  direct  S-M  interaction,  the  adsorbents  possess  high  selectivity  for 
sulfur  compounds,  but  the  steric  hindrance  make  them  difficult  to  remove  sulfur  from  DMDBT 
(dimethyl  dibenzothiophene)  etc.  In  this  project,  the  investigators  believe  that  both  of  the 
positive  effects  of  the  above  two  mechanisms  may  be  utilized  by  selecting  a  proper  adsorbent. 
Thiophene  has  two  lone  pairs  of  electrons  on  the  sulfur  atom:  one  pair  lies  on  the  six-electron  % 
system  and  the  other  lies  in  the  plane  of  the  ring.  Thiophene  can  act  either  as  an  n-type  donor  by 
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donating  the  lone  pairs  of  electrons  that  lie  in  the  plane  of  the  ring  to  the  adsorbent  (direct  S-M  a 
bond)  or  as  a  7r-type  donor  by  utilizing  the  delocalized  n  electrons  of  the  aromatic  ring  (ji  bond) 
to  form  a  7r-type  complex  with  the  metal  ions  [15].  On  the  other  hand,  the  sulfur-adsorption 
capacity  and  selectivity  of  adsorbents  can  be  further  improved  by  modifying  various  types  of 
surface  active  sites  for  sulfur  adsorption,  such  as  Lewis  acid  sites,  useful  functional  groups, 
electronic  defect  centers,  micro-structural  defects  and  so  on.  According  to  Lewis  acid-base 
theory,  most  thiophene  sulfur  compounds  in  jet  fuels  are  Lewis  base  [49],  which  are  easy  to  be 
adsorbed  at  Lewis  acid  sites.  Hence,  we  can  design  and  select  materials  that  possess  strong 
Lewis  acid  sites  to  selectively  adsorb  thiophene  sulfur  compounds  with  the  lone  pair  electrons  in 
jet  fuels.  The  Lewis  acid-base  adsorption  mechanism  is  the  interaction  between  the  acid  sites  on 
the  surface  of  adsorbents  and  thiophene  derivatives.  Additionally,  it  is  known  that  sulfur 
compounds  have  more  affinity  to  oxidation  than  their  analogue  hydrocarbons  in  jet  fuels  [35], 
therefore,  perfect  redox  properties  of  adsorbents  can  improve  oxidization  of  sulfur  compounds 
into  sulfones  and  sulfoxides.  High  conversions  of  sulfides  to  sulfones  and  sulfoxides  provide 
stronger  polarities  that  enhanced  selective  removal  of  organic  sulfur  compounds  with  solid 
adsorbents  at  ambient  temperature  and  atmospheric  pressure.  As  a  consequence,  perfect  redox 
properties  of  adsorbents  can  indirectly  increase  the  sulfur-adsorption  capacity  and  selectivity. 
Therefore,  in  the  present  study  we  developed  a  novel  hybrid  meso-porous  material  of  Ni2*  and 
Ce4+  modified  AhC^-SiC^  binary  oxide,  with  AI2O3  being  the  main  component,  for  selectively 
removing  sulfur  from  real  jet-A  fuel  with  high  sulfur  content.  The  ideal  adsorbent  model  is 
shown  in  Fig.l. 


High 
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capacity 
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Desulfurization  at 
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Easy  to  be 
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Figure  1.  Ideal  adsorbent  model  for  desulfurization  of  Jet-A  fuel 

In  the  adsorbent  of  ALCL-SiCh  with  AI2O3  being  the  main  component  oxide,  the  charge 
difference  for  each  bond  is  1/2,  and  for  all  the  bonds  the  valence  unit  is  2  according  to  Tanabe’s 
hypothesis  on  the  surface  acidity  of  binary  oxide  [50],  In  this  case,  the  Lewis  acidity  appears 
upon  the  presence  of  the  positive  charge.  Therefore,  the  big  specific  surface  and  the  Lewis 
acidity  of  the  mesoporous  Al2C>3-Si02  binary  oxide  will  greatly  enhance  adsorption  capacity  of 
thiophene  sulfurs  in  jet  fuels  at  low  temperatures.  Compared  to  conventional  adsorbent  using 
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Si02-Al2C>3  as  support  material  with  SiC>2  being  the  major  component  oxide,  the  advantage  of 
novel  adsorbent  Ni0-Ce02/Al2C>3-Si02  with  AI2O3  being  the  major  component  oxide  for 
desulfurization  at  ambient  temperature  and  pressure  has  been  reported  in  the  previous  work. 

3.2.2  Catalyst  preparation  procedures  ( information  with  patent  pending ) 

In  our  preliminary  work  [12],  adsorbent  preparation  methods — extrusion,  sol-gel,  and  wet 
impregnation,  as  shown  in  Fig.  2 — were  evaluated,  and  the  extrusion  method  was  recommended 
due  to  its  relatively  high  desulfurization  efficiency  of  the  thereby  prepared  adsorbent.  In  the 
preparation  of  the  proposed  NiO-CeCVAEC^-SiCE  adsorbent,  low  cost  raw  materials,  large  pore 
Pseudo-boehmite  (68-72%  AI2O3,  surface  area  >  300  m2/g,  pore  volume  0. 8-1.0  cm3/g)  and 
diatomite  (>95%  SiC>2),  were  respectively  used  as  the  source  of  alumina  and  silica,  and  analytical 
reagent  grade  nickel  acetate  hydrate  (Ni(CH2C00H)2.4H20)  and  cerium  acetate  hydrate 
(Ce(CH2COOH)3  5H2O)  were  employed  as  precursors.  The  mass  composition  of  Ni-Ce/ AI2O3— 
SiC>2  satisfies  that  Ni/Ce  =  10  in  mol,  Al/Si  =  15  in  mol,  and  then  Al-Si-Ox/Ni-Ce-Oy  =  9  in 
mass  (refer  to  publication  #2  under  this  project).  First,  Ni-Ce  solution  was  synthesized  by 
dissolving  nickel  acetate  hydrate  and  cerium  acetate  hydrate  in  distilled  water  under  constant 
stirring  at  60  °C.  Then  pseudo-boehmite  and  diatomite  were  mixed  with  distilled  water  and  a 
certain  amount  of  nitric  acid  in  another  container.  The  next  step  was  to  add  the  Ni-Ce  solution  to 
the  pseudo-boehmite  and  diatomite  mixture  and  stir  them  vigorously  for  10  min  to  get  the 
adsorbent  slurry  ready.  The  prepared  adsorbent  slurry  was  then  extruded  by  using  a  catalyst 
extruder  machine.  After  drying  the  extrudates  at  50  °C  in  a  forced  air  oven  overnight,  the  dried 
extrudates  were  then  calcined  in  helium  gas  in  a  tube  furnace  (MTI  OTF-1200X-80)  at  650  °C 
for  2  h.  After  cooling  down  to  the  room  temperature,  the  calcined  adsorbents  were  pelletized  and 
sieved  to  the  required  particle  size.  Before  using  for  desulfurization,  the  prepared  adsorbents 
need  to  be  dried  in  helium  gas  in  tube  furnace  at  500  °C  for  2  h  for  dehydration.  The  BET  surface 
area  and  pore  size  distribution  of  the  prepared  adsorbent  were  measured/characterized  by  a 
Micromeritics  TriStar  II  3020  equipment  using  the  N2  adsorption-desorption  method  at  liquid  N2 
temperature.  Characterization  of  the  adsorbents  helped  us  to  analyze  and  optimize  the  procedures 
and  conditions  for  fabrication  of  the  catalyst. 


Figure  2  Ni-Ce/AfOs-SiOj  adsorbents  prepared  by  three  different  methods,  (a)  Extrusion 
molding — the  best  method  adopted ;  (b)  wet  impregnation;  (c)  sol-gel. 

3.3  Desulfurization  results  from  batch  test 

Batch  tests  were  designed  to  study  the  equilibrium  and  kinetics  of  Jet-A  fuel  desulfurization 
process  based  on  the  NiO-CeC^/AECVSiCh  adsorbent.  For  each  test  run,  5  g  of  Jet-A  fuel  with  a 
measured  sulfur  concentration  of  1037  ppmw  was  put  in  a  flask  containing  pre- weighted 
adsorbent.  The  flasks  were  sealed  by  plastic  films  and  agitated  at  a  constant  rate  of  250  rpm  in  an 
orbital  shaker  (Scilogex  SK-O330-Pro)  at  room  temperature  and  atmospheric  pressure  for  24  h. 


9  |  Pag  e 


Project  final  report  August  5th,  2015 


The  desulfurization  performances  of  adsorbents  were  characterized  by  measuring  the  residual 
sulfur  concentration  in  supernatant  fuel  liquid  after  a  sufficiently  long  time  was  taken  to  reach 
the  adsorption  equilibrium.  Figure  3  shows  a  visual  comparison  of  the  commercial  Jet-A  fuel  and 
the  treated  fuel  in  our  preliminary  study.  The  color  of  the  fuel  changed  from  a  straw  color  to 
colorless,  which  proves  that  the  developed  adsorbent  has  an  impressive  selectivity  of 
organosulfur  compound.  The  total  sulfur  concentration  was  measured  by  a  Thermo  TS3000  total 
sulfur  analyzer,  which  has  a  working  range  of  0.02-5000  ppmw,  and  an  uncertainty  of  less  than 
5%  of  the  measured  value.  The  desulfurization  efficiency  (%)  and  equilibrium  adsorption 
loading  q  in  the  adsorbent  (mg/g)  were  calculated  as  follows: 

Desulfurization  efficiency  (%)  =(C0 -Cs )/Coxl00  (1) 

Loading  q  (mg-S/g-ads)  =mfue,  (C0  -Cs)/ mads  (2) 

where  Co  is  the  initial  sulfur  concentration  (ppmw),  Cs  is  the  residual  sulfur  concentration 
(ppmw),  nifuet  (g)  is  the  mass  of  fuel  sample,  and  triads  (g)  is  the  mass  of  the  adsorbent,  which 
includes  the  catalyst  and  support  material. 


Figure  3  Left:  original  Jet-A  fuel;  Right:  treated  Jet-A  fuel 


3.3.1  Effects  of  adsorbent  calcination  temperature  and  time 

To  prepare  the  adsorbent,  any  dried  adsorbent  extrudate  was  calcined  at  a  constant 
temperature  for  3  hours  in  a  helium  gas  atmosphere.  In  our  preliminary  study  helium  was  proved 
to  be  a  better  calcination  gas  than  nitrogen,  argon  and  air.  Different  temperatures  ranging  from 
400-  800  °C  were  chosen  for  the  calcination  process,  and  the  prepared  adsorbent  was  tested  for 
desulfurization  efficiency  in  order  to  optimize  the  calcination  temperature.  For  each  test, 
sufficient  mass  of  4  g  adsorbent  was  put  into  5g  Jet-A  fuel,  and  the  sulfur  concentration  in  the 
Jet-A  fuel  sample  was  measured  after  24  h  adsorption. 

The  desulfurization  efficiencies  of  each  adsorbent  calcined  at  different  temperatures  are 
shown  in  Fig.  4.  It  is  evident  that  the  desulfurization  performance  of  the  adsorbent  significantly 
improves  as  the  calcination  temperature  increases  from  400  °C  to  650  °C.  From  650  °C  to  800  °C 
the  desulfurization  performance  slightly  decreases  with  the  increase  of  calcination  temperature. 
The  increase  of  calcination  temperature  accelerates  adsorbent  sintering,  promoting  the 
combustion  of  organic  matters  and  the  volatilization  of  physically  and  chemically  combined 
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water,  which  enhances  the  surface  activation  energy  of  the  adsorbents.  However,  extremely  high 
temperatures  can  induce  over-sintered  adsorbents,  thus  lowering  the  adsorption  activity.  In 
conclusion,  the  adsorbent  calcined  at  650  °C  has  the  best  desulfurization  efficiency  of  96.5%  on 
Jet-A  fuel  with  a  1,037  ppmw  initial  sulfur  concentration. 


Figure  4  Effect  of  calcination  temperature  on  adsorbent  desuljurization  performance 

Besides  the  calcination  temperature,  the  calcination  time  can  also  influence  the  adsorbent 
desulfurization  performance.  Figure  5  presents  the  effect  of  calcination  time  on  adsorbent 
desulfurization  efficiency  for  Jet-A  fuel.  All  adsorbents  were  calcined  at  650  °C  in  a  helium  gas 
atmosphere.  Again,  for  each  of  the  tests,  sufficient  adsorbent  and  time  were  used  and  sulfur 
concentration  could  no  longer  be  reduced  by  adding  more  adsorbent.  It  shows  that  the  variation 
of  calcination  time  in  0.5  to  4  hours  affects  the  desulfurization  efficiency  about  3%,  which  is  less 
significant  compared  to  the  effect  of  the  calcination  temperatures.  Nevertheless,  the  experimental 
results  showed  that  2  hours  is  the  optimal  calcination  time,  which  is  used  to  prepare  the 
adsorbent  samples  in  our  study.  The  maximum  desulfurization  efficiency  could  reach  up  to  98%. 


Figure  5  Effect  of  calcination  time  on  adsorbent ’s  desulfurization  performance 

3.3.2  Preliminary  adsorbent  screening 

The  adsorbents  with  different  proportions  of  metal  catalysts  and  support  materials  were 
prepared  under  the  same  procedure,  followed  by  calcination  at  650  °C  for  2  h  in  helium  gas. 
There  were  10  samples  prepared  in  this  study.  Their  molar  ratios  of  A1  to  Si  for  the  support,  Ni  to 
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Ce  for  the  catalyst,  and  mass  ratio  of  Al-Si-Ox  to  Ni-Ce-Oy,  varied  in  the  ranges  of  10-20,  5-15, 
and  4-11,  respectively.  For  each  test,  sufficient  mass  of  4  g  adsorbent  was  put  into  5g  Jet-A  fuel, 
and  the  sulfur  concentration  in  the  Jet-A  fuel  sample  was  measured  after  24  h  adsorption.  The 

desulfurization  efficiencies  of  the  samples  are  compared  in  Figure  6.  Sufficient  adsorbent  and 
time  was  used  and  sulfur  concentration  could  no  longer  be  reduced  by  adding  more  adsorbent. 
The  curves  in  Fig.  6  indicate  that  the  performance  of  the  10  samples  varies  significantly.  It  is 
understood  that  the  disparate  loading  of  metal  catalysts  on  support  materials  results  in  different 
specific  surface  areas,  as  well  as  different  pore  diameters  and  volumes  on  the  adsorbent,  which 
drastically  affects  the  sulfur  adsorption  [51,52],  The  highest  desulfurization  efficiency  reaches 
98%.  The  composition  of  the  metal  catalyst  and  support  materials  from  this  sample  with  the  best 
performance  was  used  in  the  following  further  study.  The  adsorbent  sample  with  the  best 
desulfurization  performance  had  Ni/Ce=10  in  mole,  Al/Si=15  in  mole,  and  Ni-Ce-Ox/Al-Si- 
Oy=9  in  mass. 


Adsorbent  sample  number 

Figure  6  Effect  of  composition  on  adsorbent  desulfurization  performance 

3.3.3  Effect  of  adsorbent-to-fuel  mass  ratio 

The  effect  of  the  adsorbent  mass  in  a  certain  amount  of  fuel  was  investigated  by  varying  the 

adsorbent  from  0  to  8  g  for  5  g  Jet-A  fuel  in  the  desulfurization  tests.  After  shaking  the  test 
samples  with  the  fuel  on  a  vibration  plat  with  250  rpm  for  24  h,  the  sulfur  concentrations  in  the 
fuel  were  measured,  and  the  results  are  shown  in  Fig.  7.  Obviously,  with  the  increase  of  the  mass 
of  adsorbent  from  0  to  4  g,  less  and  less  sulfur  compounds  were  left  in  the  fuel,  which  is 
reflected  as  a  significant  increase  of  the  desulfurization  efficiency.  When  the  mass  of  the 
adsorbent  is  more  than  4  g,  the  sulfur  concentrations  of  the  after-treated  fuel  does  not  decrease 
any  more,  w'hich  means  there  is  no  increase  of  the  desulfurization  efficiency. 
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Figure  7  Effect  of  adsorbent  mass  on  adsorbent  desulfurization  performance  (at  room 

temperature ) 

The  increase  of  sulfur  adsorption  shown  in  the  curve  is  due  to  the  increase  of  the  mass  of  the 
absorbent  from  0  to  4  g,  which  actually  increases  the  active  surface  area,  or  sites  for  sulfur 
adsorption.  When  the  mass  of  the  adsorbent  is  more  than  4  g,  the  adsorption  equilibrium  between 
adsorbate  and  adsorbent  is  reached  so  that  the  desulfurization  efficiency  is  almost  constant.  This 
effect  is  the  so  called  overcrowding  of  particles  as  elaborated  in  the  reference  [9]. 

3.3.4  Kinetic  models  and  equations 

In  order  to  better  understand  the  adsorption  process,  the  commonly  used  kinetic  models 
including  the  pseudo-first  order  model,  the  pseudo-second  order  model,  and  the  intraparticle 
diffusion  model  were  examined  to  best  describe  the  sulfur  adsorption  kinetics  of  the  new 
developed  adsorbents. 

The  pseudo-first  order  Lagergren  model  is  expressed  as 

rq=kf{qe-qt)  =  dqtldt  (3) 

where  qt  is  the  load  or  uptake  of  sulfur  per  unit  mass  of  adsorbent  in  a  time  period  of  t,  rq  is  the 
sulfur  adsorption  rate  in  mg. gamin'1,  kf  is  the  pseudo-first  order  model  rate  constant  (min'1),  and 
qe  is  the  equilibrium  sorption  uptake  (mg.g'1).  From  kinetic  batch  test  of  sulfur  loading  versus 
adsorption  time,  the  experimental  equilibrium  sorption  uptake  qe  was  determined  at  the 
adsorption  time  that  is long  enough  to  reachtheequilibriumstate.  And  the  calculated  value  of  qe_ 

can  be  used  to  compare  with  the  experimental  value  to  validate  the  kinetic  model.  In  this  test  a 

sufficient  mass  of  4  g  adsorbent  was  put  into  5  g  Jet-A  fuel  to  study  the  sulfur  adsorption  kinetics. 

The  experimental  equilibrium  sorption  uptake  ge  was  1.321  mg/g  at  the  adsorption  time  of  22  h , 
The  integrated  form  of  Eq.  (3)  is  given  as 


Hqe-q,)=^qe~k/t  (4) 

The  pseudo-first  order  model  inherently  assumes  the  sulfur  concentration  in  the  ambient  (which 

refers  to  the  5  g  Jet-A  fuel  in  this  test)  doesn’t  change  with  time,  and  l/kfis  a  time  constant  in  the 
above  equation.  However,  in  this  test  the  sulfur  concentration  Cs  in  the  Jet-A  fuel  sample 

reduced  significantly  as  the  adsorption  time  increased,  because  the  mass  of  the  Jet-A  fuel  was 

limited  and  sufficient  mass  of  adsorbent  was  added  into  the  fuel.  At  time  f,  the  mass  balance  of 

sulfur  is  expressed  as, 

mf(C0-Cs)  =  m,,0  ~  mfCs  =  q,  ■  Wads  (5) 

where  Co  is  the  initial  sulfur  concentration  (ppmw),  Cs  is  the  residual  sulfur  concentration 

(ppmw),  mjuej  (g)  is  the  mass  of  fuel  sample,  m,^  (g)  is  the  mass  ofthe  adsorbent,  and  m^p  (g)  is 

the  total  mass  of  the  sulfur.  At  equilibrium  state,  there  was  still  sulfur  left  in  the  fuel  because  of 

the  balance  of  adsorption  and  desorption. 

ms,Q=mfCe+qemads  (6) 

where  Ct  is  the  sulfur  concentration  in  the  Jet-A  fuel  at  equilibrium  state.  By  substituting  Eq.  (6) 

to  Eq.  (5),  the  following  equation  is  obtained  as: 
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mf  , 

— -  (C,-C.) 


m 


ads 


The  pseudo-second  order  model  is  obtained  by  combining  Eq.  (3)  and  Eq.  (7): 

2 

=dq,/dt 


where  ks  is  the  pseudo-second  order  model  rate  constant  (g-mg^min-1), 


K^— 


m  f  cs-ce 


The  integrated  form  of  Eq.  (8)  is: 


q,  ksq2e  qe 


(7) 


(8) 

(9) 

(10) 


The  intraparticle  diffusion  model  is  presented  in  the  form  of 

q,=k,tm+C  (11) 

where  kj  is  the  intraparticle  diffusion  model  rate  constant  (mg.g’lmin"'lS)  and  C  is  a  constant 
(mg.g'1). 

From  the  experimentally  measured  data  of  qt  versus  time  t,  linear  regression  of  Eqs.  (4),  (10) 
and  (11)  can  be  conducted,  and  thus  the  qe,  k  n  ks  and  C  in  the  equations  have  been  found  as 

given  in  Table  2. 


_ Table  2  Kinetic^  analysis  results^  ofthree  difiererU  models _ _ 

Pseudo-first  order  model  Correlation  coefficient  R‘ 

qe  cai  =0.1078  (mg/g)  kf  =0.00232  (min'1)  0.9584 

Pseudo-second  order  model 

qe,cai  =1.321  (mg/g)  ks  =0.7569  (g.mg^min1)  0.9999 

Introparticle  diffusion  model 

C=  1.1467  (mg/g) k,  =0.0058  (mg.g 'min' °'5) 0.5851 

The  kinetic  analysis  results  following  the  above-mentioned  three  different  models  are  shown 
in  Fig.  8  to  Fig.  10.  It  was  found  that  the  pseudo-second  order  model  can  best  fit  the  current 
experimental  data.  Table  2  lists  the  data-fitted  qe  and  the  correlation  coefficient  R2  for  all  three 
models.  Because  of  the  excellent  fit  of  the  experimental  data  to  the  equation  from  the  pseudo- 
second  order  model,  the  thus  obtained  equilibrium  loading  qe  agrees  with  the  experimental  data 
in  Fig.  7  very  well.  This  leads  to  the  conclusion  that  the  kinetics  of  the  desulfurization  of  Jet-A 
fuel  using  the  currently  developed  adsorbent  can  be  well  described  by  the  pseudo-second  order 
model. 

The  intraparticle  diffusion  model  has  a  poor  fitting  with  the  experimental  data,  which 
indicates  that  the  diffusion  process  may  have  more  than  one  rate-controlled  stage,  which  needs 
two  piecewise  data  regressions.  This  is  clearly  illustrated  by  plotting  the  adsorbent  loading 
versus  t  as  shown  in  Fig.  8.  Kumar  et  al.  [9]  also  proposed  that  the  adsorption  rate  is  controlled 
by  both  chemical  adsorption  and  intraparticle  transport  at  the  beginning,  and  then  followed  by 
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meso-pore  diffusion  and  micro-pore  diffusion. 


Figure  8  Kinetic  modeling  analysis  of  pseudo-first  order  model 


Figure  9  Kinetic  modeling  analysis  of  pseudo-second  order  model 


Figure  10  Kinetic  modeling  analysis  of  intraparticle  diffusion  model 
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Figure  11  Plot  of  loading  versus  t 


3.3.5  Equilibrium  isotherm 

The  equilibrium  adsorption  loads  of  sulfur  on  the  adsorbent  from  the  present  experimental 
data  are  shown  in  Fig.  12  in  respect  of  the  loading,  q,  against  the  residual  sulfur  concentration, 
Cs.  The  curves  are  also  called  isotherm.  Two-parameter  Freundlich  and  Langmuir  isotherm 
equations  and  three-parameter  Redlich-Peterson  and  Brunauer-Emmett-Teller  (BET)  isotherm 
equations  were  used  to  analyze  the  adsorption  data.  The  linear  forms  of  the  four  equations  are 
given  as 


Freundlich: 

Langmuir: 

Redlich-Peterson: 

BET: 


In  q  =  In  KF  +  —  In  Cs 
n 

1  __L  1  1 


d  dm  K,qm  Cs 


In 


(  „  a 


kr-~  1 


=  (3\nCs  +  \naR 


C 


Ka-\C .  1 


 -“-a 


+  - 


q{Cm~Cs)  qmKBCm  qmK 


(12) 

(13) 


(14) 


(15) 


The  regression  of  the  experimental  data  in  Fig.  12  following  Eq.  (12)  allowed  us  to  find  Kf 
and  1/n.,  and  similarly  we  can  find  K  L  and  qm  for  Eq.  (13).  The  data  regression  for  Eq.  (14) 
requires  us  to  guess  a  KR  and  find  /?  and  aR  for  the  best  fit  of  the  data.  Through  trail-and-error 
analysis,  the  current  analysis  found  the  best  fit  at  a  correlation  coefficient  R2  =0.9789.  The  data 
regression  for  Eq.  (15)  requires  us  to  guess  a  value  for  Cm  and  then  find  qm  and  KB  for  the  best 
fit  to  the  data.  Through  trail-and-error  analysis,  the  current  analysis  found  the  best  fit  at  a 
correlation  coefficient  R2  =0.9986. 
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In  order  to  better  estimate  the  deviation  of  the  calculated  values  from  data-fitted  equations 
and  experimental  data,  the  sum  of  the  square  of  errors  (SSE)  is  calculated, 

SSE  =  YX%cal-<l,,exP)  (16) 

i 

The  Marquardt’s  percent  standard  deviation  (MPSD)  error  function  is  also  calculated  from  the 
equation 


MPSD  =  100, 


1  Y, 

r  _  \2 

7 1  xal  h  l.exp 

nm  —nni 

m  p  1 

^  Qi,exp  j 

(17) 


where  the  subscript  ‘exp’  and  ‘cal’  represent  the  experimental  data  and  calculated  value  from  the 
data-fitted  equation,  respectively,  nm  is  the  total  number  of  measurement  points  and  np  is  the 
number  of  to-be-determined  parameters  in  the  model  equations. 

The  curves  from  the  data-fitted  equations  for  the  Freundlich,  Langmuir,  Redlich-Peterson, 
and  BET  models  are  also  shown  in  Fig.  12.  The  obtained  parameters  from  data  regression,  the 
correlation  coefficient  R2 ,  SSE,  and  MPSD  values  of  each  model  are  given  in  Table  3.  It  is  seen 
that  all  four  models  can  more  or  less  fit  the  experimental  data  in  the  investigated  range  of 
residual  concentration.  A  very  good  fit  is  seen  for  the  BET  model,  which  has  the  highest 
correlation  coefficient  R2,  and  the  lowest  SSE  and  MPSD,  which  indicates  that  the  BET  model  is 
the  most  suitable  one.  The  BET  model  was  believed  most  suitable  for  the  case  of  multi-layer 
adsorption  and  thus  high  desulfurization  efficiency.  The  good  fit  of  the  BET  model  to  the  current 
data  indicates  that  the  adsorbent  developed  offers  multi-layer  adsorption.  The  satisfactory  fitting 
of  the  equations  with  the  experimental  data  also  indicates  that  the  organosulfur  compounds  in 
Jet-A  fuel  can  be  well  represented  by  the  total  sulfur  amount  in  the  fuel. 


Sulfur  concentration  (ppmw) 

Figure  12  Adsorption  equilibrium  isotherm  at  room  temperature 


Table  3  Isotherm  parametersjor  sulfur  adsorption 


Kf 

1/n 

Freundlich  Equation 

R2 

SSE 

MPSD 

[(mg/g)/(mg/kg)1/n] 

0.4007 

0.9513 

0.5001 

16.747 

0.2998 

Kl  (kg/mg) 

qm  (mg/g) 

Langmuir  Equation 

R2 

SSE 

MPSD 

0.0158 

3.9736 

0.9609 

0.2897 

11.682 

Redlich-Peterson  Equation 
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Kr  (kg/g) 

aR  (kg/mg) 

P 

R2 

SSE 

MPSD 

0.219 

0.4684 

0.6657 

BET  Equation 

0.9789 

0.3074 

14.985 

Cm  (mg/kg) 

qm  (mg/g) 

Kb 

R2 

SSE 

MPSD 

6250 

3.6426 

103.4261 

0.9986 

0.1435 

11.5193 

3.4  Desulfurization  results  from  fixed  bed  reactor 

Fixed  bed  tests  were  conducted  for  dynamic  sulfur  adsorption  measurements  and  investigated 
effects  of  fuel  flow  rate,  adsorbent  particle  size,  fixed  bed  dimensions,  etc.  [48].  A  vertically- 
oriented  stainless  steel  column  fully  packed  with  adsorbents  was  used  as  the  fixed  bed  reactor. 
The  experimental  system  consisted  of  a  HPLC  pump  (high  performance  liquid  chromatography 
pump),  a  fuel  tank,  and  the  fixed  bed  column.  Nitrogen  was  used  to  flush  out  the  residual  fuel 
after  each  test.  Figure  13  shows  schematically  the  test  rig  of  the  fixed  bed  adsorption  system. 
After  the  tubing  was  rinsed  with  distilled  water  and  dried,  pre-weighted  adsorbent  was  packed 
into  the  fixed  bed  column.  The  HPLC  pump  (SSI  B2300)  was  used  to  feed  the  Jet-A  fuel  to  the 
adsorption  bed  at  the  desired  flow  rate.  The  effluent  fuel  from  the  adsorption  bed  was  sampled 
every'  15  min.  and  the  sulfur  concentration,  C(t),  of  every  sample  at  a  time  was  measured  by  the 
Thermo  TS  3000  total  sulfur  analyzer.  The  treated  accumulated  fuel  after  desulfurization  was 
collected  in  a  beaker.  The  sulfur  adsorption  capacity  of  the  adsorbent  was  calculated  as  follows: 

Sulfur  adsorption  capacity  q= — — — - — f  (Cn-C(t))dt  (18) 

c  1000-mfl£/i  J°v  v  ,r 

where  qc  has  a  unit  of  mg-S/g-adsorbent,  p/(g/ml)  is  the  density  of  fuel,  V/  (ml/min)  is  the  fuel 
flow  rate,  ma(js  (g)  is  the  mass  of  packed  adsorbent,  and  Co  (ppmw)  is  the  initial  sulfur 
concentration  in  the  Jet-A  fuel. 


Figure  13  Schematic  of  the  fixed-bed  sulfur  adsorption  system 

All  the  results  were  obtained  with  repeatability  validation.  We  repeated  and  prepared  the 
absorbent  in  two  times.  Also,  each  test  was  run  twice  in  the  data  collection  stage  in  order  to 
confirm  the  repeatability  of  experiments. 

3.4.1  Effect  of  LHSV  (liquid  hourly  space  velocity) 

The  desulfurization  performance  of  adsorbent  strongly  depends  on  the  LHSV  (liquid  hourly 
space  velocity),  which  can  be  changed  by  using  different  flow  rate  of  fuel.  In  order  to 
demonstrate  this  effect  for  the  current  adsorbent,  adsorptive  desulfurization  experiments  were 
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carried  out  under  flow  rates  of  0.1,  0.2,  0.3,  and  0.4  ml/min  ( in  corresponding  LHSVs  of  0.168, 
0.336,  0.504.  and  0.672  h  ,  respectively  )  while  keeping  other  parameters  unchanged. 
Adsorbents  with  particle  size  between  0.5  and  3  mm  were  used  in  each  test.  The  adsorption 
capacities  versus  the  mean  sulfur  concentration  of  the  treated  (accumulated)  fuel  at  different  flow 
rates  are  shown  in  Fig.  14.  It  is  seen  that  as  the  flow  rate  decreases  from  0.4  to  0.1  ml/min,  the 
adsorbent  capacity  increases  significantly.  When  the  accumulated  treated  fuel  sulfur 
concentration  reached  30  ppmw,  the  capacity  of  sulfur  adsorption  at  the  flow  rate  of  0.1  ml/min 
is  about  two  times  of  that  of  the  0.4  ml/min  flow  rate.  This  is  because  that  at  the  high  flow  rate, 
the  per  unit  volume  of  fuel  contacts  the  adsorbent  bed  in  a  shorter  time,  therefore  the  sulfur 
compounds  do  not  have  sufficienct  time  to  diffuse  into  the  pores  of  the  adsorbent.  After  all, 
although  setting  a  low  flow  rate  increases  the  desulfurization  performance  of  the  adsorbent,  the 
fuel  desulfurization  rate  can  be  too  small,  particularly  when  the  flow  rate  is  less  than  0.1ml  /min. 
Therefore,  in  the  tests  hereafter,  the  fuel  flow  rate  was  selected  to  be  no  less  than  0.1  ml/min. 


Figure  14  Effect  of  flow  rate  on  sulfur  removal  performance 

3.4.2  Effect  of  adsorbent  particle  size 

To  study  the  effect  of  particle  size  on  the  sulfur  adsorption  capacity,  three  experiments  were 
carried  out  under  the  same  conditions  except  for  the  particle  size.  The  test  results  are  shown  in 
Figure  15.  Since  the  total  mass  of  the  adsorbent  is  the  same,  the  packed  bed  has  different 
porosities  due  to  the  difference  of  the  particle  sizes.  Figure  15  shows  that  with  smaller  particle 
size  a  higher  sulfur  adsorption  capacity  was  obtained.  Because  of  this,  the  lifetime  of  adsorbent 
is  also  significantly  increased.  This  is  understandable  because  with  a  smaller  particle  size  for  the 
adsorbent  the  sulfur  mass  diffusion  can  reach  the  center  of  the  particles  easier  and  thus  increase 
the  utilization  of  the  adsorbent  material. 


Sulfur  content  in  effluent  fuel  (ppmw) 

Figure  15  Effect  of  adsorbent  particle  size  on  sulfur  removal  performance 
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3.4.3  Effect  of  fixed  bed  dimensions 

While  keeping  the  same  volume  of  adsorbent,  the  tube  diameter  and  the  length  of  the  packed 
bed  can  be  varied.  To  investigate  the  effect  of  the  packed  bed  dimensions  on  the  desulfurization 
capacity  of  the  adsorbent,  three  different  diameters  and  lengths  of  the  packed  bed  were  chosen 
and  filled  with  adsorbent  of  particle  sizes  0-125  pm.  The  fuel  feeding  flow  rate  and  the  mass  of 
packed  adsorbent  were  identical  in  these  tests.  Figure  16  shows  the  results  of  sulfur  adsorption 
capacity  versus  the  total  sulfur  concentration  of  the  accumulated  fuel  after  desulfurization. 

The  test  results  show  an  optimal  diameter  for  the  fixed  volume  of  the  reactor.  Of  the  tested 
three  diameters,  the  packed  bed  with  a  diameter  of  7.75  mm  gave  the  highest  adsorption  capacity 
of  1.55  mg  S/g  adsorbent  at  50  ppmw  S  in  effluent  fuel.  When  the  volume  of  the  adsorbent  bed 
and  the  fuel  flow  rate  were  both  unchanged  in  the  tests,  the  difference  of  the  packed  bed  diameter 
only  causes  different  velocities  of  the  fuel  when  it  flows  through  the  packed  bed.  On  the  other 
hand,  the  time  of  the  fuel  passing  through  the  packed  bed  is  the  same  for  all  the  cases.  Therefore, 
the  different  adsorption  capacity  in  the  three  cases  is  only  due  to  the  difference  of  flow  velocities. 
A  very  low  velocity  does  not  contribute  to  a  relatively  good  mass  diffusion/adsorption  and, 
therefore,  we  see  that  the  highest  sulfur  adsorption  capacity  is  from  the  case  of  a  packed  bed 
diameter  of  7.75  mm. 


Sulfur  content  In  effluent  fuel  (ppmw) 

Figure  16  Effect  of  fixed  bed  dimensions  on  sulfur  removal  performance 

A  complete  adsorption  curve,  effluent  sulfur  concentration  versus  adsorption  capacity  is 
shown  in  Fig.  17,  which  was  obtained  at  conditions  of  L=203.2  mm,  D=7.747  mm,  fuel  flow 
rate=0.1  ml/min,  LHSV=0.63,  adsorbent  particle  size=0~0.125  mm,  fuel  initial  sulfur 
concentration^, 037  ppmw,  mass  of  adsorbent=4.88  g.  At  a  breakthrough  sulfur  concentration  of 
10  ppmw  a  very  high  sulfur  adsorption  capacity  of  0.633  mg  S/g  adsorbent  was  achieved,  which 
is  significantly  higher  than  the  value  of  0.12  mg  S/g  adsorbent  reported  in  literature  [53].  The 
sulfur  adsorption  capacity  reached  1 .98  mg  S/g  adsorbent  at  the  point  of  30  ppmw  of  mean  sulfur 
concentration  in  the  treated  accumulated  fuel. 
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Capacity  (mg  S/g  adsorbent) 

Figure  17  A  complete  adsorption  curve 

(L=203.2  mm;  D=7. 747  mm:  fuel  flow  rate=0. 1  ml/min;  LHSV=0.63;  adsorbent  particle 
size=0~0.1 25  mm;  fuel  initial  sulfur  concentration= 1 037 ppmw;  adsorbent  mass=4.88  g) 

3.4.3  Effect  of  the  amount  of  nitric  acid  used  for  absorbent 

In  the  adsorbent  preparation  process,  a  small  amount  of  nitric  acid  was  added  in  order  to  help 
the  alumina  sol  formation  when  preparing  the  adsorbent  support.  However,  it  w'as  observed  that 
the  amount  of  nitric  acid  had  some  effect  in  the  adsorbent  desulfurization  performance.  In  this 
study,  the  total  128.25  gram  of  prepared  mixture  of  pseudo-boehmite,  diatomite,  cerium  acetate 
hydrate,  nickel  acetate  hydrate,  and  distilled  water,  was  added  with  different  amount  of  nitric 
acid  from  0  to  5  g,  and  the  desulfurization  performance  of  the  processed  adsorbent  was  tested. 
Sulfur  adsorption  capacities  of  the  five  adsorbents  with  0  g,  1.5  g,  3  g,  4  g  and  5  g  nitric  acid 
versus  total  sulfur  concentration  in  effluent  fuel  are  shown  in  Fig.  18.  Breakthrough  points  for  all 
five  adsorbents  occurred  at  around  50  ppmw  S  in  effluent  fuel,  and  it  indicated  that  the  less  the 
nitric  acid  was  added  in  the  adsorbent  preparation  stage,  the  better  desulfurization  performance 
the  adsorbent  had.  At  50  ppmw  S  breakthrough  point,  the  sulfur  adsorption  capacity  for 
adsorbent  with  no  use  of  nitric  acid  was  2.95  mg-S/g-ads,  while  the  obtained  adsorption  capacity 
was  2.5  mg-S/g-ads  for  adsorbent  with  1.5  g  nitric  acid,  2.2  mg-S/g-ads  for  adsorbent  with  3  g 
nitric  acid,  2.16  mg-S/g-ads  for  adsorbent  with  4  g  nitric  acid  and  2.1  mg-S/g-ads  for  adsorbent 
with  5  g  nitric  acid.  Although  the  difference  of  capacities  for  the  five  adsorbents  was  not 
significant  at  50  ppmw  S  in  the  effluent  fuel,  the  gap  became  larger  as  the  sulfur  concentration  in 
the  effluent  fuel  increases. 


Figure  18  Effect  of  acid  on  sulfur  removal  performance 
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Whereas  the  mechanism  of  desulfurization  by  using  the  selected  metal  and  support  materials 
from  the  chemistry  point  of  view  has  been  discussed,  some  physical  details  of  the  adsorbent 
should  be  interesting  to  examine.  The  BET  surface  areas  were  256,  253,  252,  249  and  249  m2/g 
for  adsorbents  containing  0,  1.5,  3,  4  and  5  g  nitric  acid,  respectively.  The  pore  size  distributions 
of  the  five  adsorbents  are  shown  in  Fig.  19.  The  surface  areas  slightly  decrease  with  the 
increased  addition  of  nitric  acid.  Also,  the  addition  of  nitric  acid  in  adsorbent  preparation  reduces 
the  pore  size  of  the  adsorbent.  The  BJH  average  pore  diameter  (4V/A)  were  92,  74,  72,  70  and  65 
A  for  adsorbents  with  0  g,  1.5  g,  3  g,  4  g  and  5  g  nitric  acid,  respectively.  The  result  of  reduction 
of  both  surface  area  and  pore  size  is  correlated  to  the  reduction  of  desulfurization  performance  of 
the  adsorbent.  It  is  understood  that  the  increase  of  surface  area  of  the  adsorbent  has  positive 
effect  to  the  amount  of  adsorption  of  sulfur  compound,  whereas  the  large  pore  size  in  the 
adsorbent  is  favorable  for  the  mass  transfer  resistance  in  the  adsorbent,  which  can  be  even  more 
important.  Sulfur  compound  with  large  molecules  may  not  able  to  pass  through  the  small  pores 
and  thus  the  adsorption  efficiency  can  be  low.  The  reduction  of  BET  surface  area  and  BJH 
average  pore  diameter  caused  by  addition  of  acid  might  be  due  to  the  fact  that  nitric  acid  could 
dissolve  pseudo-boehmite  to  prepare  alumina  sol.  More  pseudo-boehmite  was  dissolved  as  more 
nitric  acid  was  added  in  the  adsorbent  preparation.  Because  the  structure  of  alumnia  sol  is 
different  to  that  of  pseudo-boehmite,  the  addition  of  nitric  acid  caused  structure  change  of  the 
prepared  adsorbents  and  so  as  the  desulfurization  performance  of  the  adsorbents  in  fixed  bed 
tests. 


Logarithm  of  pore  width  W  (A) 

Figure  19  Pore  size  distribution  of  adsorbents  with  various  acid  amounts 

(Unit  of  pore  size  W:  A) 


3.4.4  Effect  of  compositions  of  materials  in  adsorbent 

The  loading  ratio  of  the  two  metal  oxides  and  the  ratio  of  alumina  to  silica  affect  the 
dispersion  of  active  sites  and  interactions  of  metals  with  support  materials.  As  the  result,  the 
ratios  of  Ni/Ce  and  Al/Si  could  significantly  influence  the  desulfurization  performance. 
Obviously,  the  desulfurization  mechanism  includes  both  chemical  and  physical  aspects.  The 
following  measured  physical  properties  of  the  adsorbent  will  help  for  better  understanding  of  the 
desulfurization  performance  of  different  adsorbent  samples  from  mass  transfer  and  surface 
adsorption  point  of  view. 


Table  4  Molar  ratios^  in  meted  oxide  materials  and  properties  of  adsorbent 


Adsorbent 


Ni/Ce 
in  mole 


Al/Si 

in 


Al-Si-Ox/Ni- 

Ce-Oy 


Total  wt.% 
of 


Surface 

area 


Average  pore 
size 
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mole 

in  mass 

metal  oxide 

(m  /g) 

(A) 

Support 

n/a 

15 

9 

0 

288 

105 

Ads-1 

10 

15 

10 

6.6 

277 

93 

Ads-2 

10 

15 

9 

7.4 

256 

92 

Ads-3 

8 

15 

9 

7.3 

242 

96 

Ads-4 

12 

15 

9 

8.2 

235 

91 

Ads-5 

10 

15 

8 

8.0 

266 

86 

Ads-6 

10 

12 

9 

7.4 

245 

94 

In  order  to  observe  the  optimized  adsorbent  performance  in  fixed  bed  reactor  test,  six 
adsorbents  labeled  as  Ads-1  to  Ads-6  were  studied  under  the  same  fixed  bed  operation  condition. 
And  the  performance  of  pure  support  material  was  also  tested  as  a  baseline  sample.  For  Ads-1  to 
Ads-5,  the  support  materials  were  the  same,  meaning  the  ratio  of  pseudo-boehmite  versus 
diatomite  are  the  same;  whereas  the  loading  of  NiO  varies  from  5  wt.%  to  7  wt.%,  and  the 
loading  of  CeC>2  varies  from  1.2  wt.%  to  1.6  wt.%.  The  total  metal  oxide  loading  is  given  in 
Table  4.  The  samples  Ads-6  and  Ads-2  had  the  same  metal  oxide  loadings  for  both  the  ratios  of 
Ni0/Ce02  and  metal  materials  versus  support  materials.  But,  in  the  support  material  of  Ads-6 
the  ratio  of  pseudo-boehmite  versus  diatomite  was  different  from  that  of  sample  of  Ads-2.  Table 
4  also  gives  the  mole  ratios  of  Ni  versus  Ce,  the  total  weight  percentage  of  metal  oxide  in  the 
adsorbent,  as  well  as  the  measured  surface  area  and  pore  size. 

The  BET  measurement  showed  that  the  surface  area  of  Ads-1  to  Ads-4  decreases  from  277 
m2/g  to  235  m2/g  as  the  metal  oxide  loading  increases.  If  only  consider  support  materials  with  no 
metal  oxide  loading,  the  surface  area  was  measured  to  be  288  m2/g.  This  means  that  the  loading 
of  metal  oxide  actually  slightly  reduces  the  surface  area  of  the  adsorbent.  The  pore  size 
distribution  curves  of  Ads-1  to  Ads-6  didn't  vary  significantly  as  shown  in  Fig.  20.  The  BJH 
average  pore  diameter  for  only  support  material  and  Ads-1  to  Ads-6  were  105,  93,  92,  96,  91,  86 
and  94  A,  respectively.  The  diffusion  limitation  of  large  sulfur  compound  molecule  was  the  key 
factor  in  the  sulfur  adsorption,  and  adsorbent  with  large  pores  can  have  minimized  diffusion 
limitation  and  thus  increase  sulfur  adsorption  performance.  The  pore  size  increase  is  more 
effective  on  adsorbent  sulfur  adsorption  capacity  than  the  BET  surface  area  increase.  For  the 
currently  developed  NiO-CeCVAEC^-SiC^  adsorbent,  the  measured  pores  were  mainly 
mesopores  (pore  width  between  2-50  nm),  and  the  BJH  average  pore  diameter  4V/A  (~90  A)  of 
the  new  adsorbents  are  much  larger  than  the  BJH  average  pore  size  of  adsorbents  made  of  Y- 
zeolite  (~20  A)  [14],  MCM-41  (~30  A),  and  SBA-15  (~60  A). 
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Logarithum  of  pore  width  W  (A) 

Figure  20  Pore  size  distribution  of  Ads-1  to  Ads-6  and  pure  support 
(Unit  of  pore  size  W:  A) 


Sulfur  content  in  effluent  fuel  (ppmw) 

Figure  21  Capacity  versus  sulfur  content  in  effluent  fuel  for  various  adsorbents 
(Co:  1140  ppmw;  LHSV:  0.63  h'1 ;  particle  size:  0-125  pm) 


The  pore  size,  surface  area,  and  the  total  weight  percentage  of  metal  oxide  in  the  adsorbent 
work  together  to  influence  the  desulfurization  performance.  The  adsorption  capacities  of  Ads-1 
to  Ads-6  against  the  sulfur  content  in  effluent  fuel  are  shown  in  Figure  21.  While  larger  surface 
area  is  favorable;  the  total  weight  concentration  of  metal  oxide  and  the  pore  size  also  play 
significant  roles.  Overall,  Ads-2  had  the  best  desulfurization  performance  than  other  adsorbents. 
Its  sulfur  adsorption  capacities  were  1.22  and  2.95  mg-S/g-ads  at  10  and  50  ppmw  S  in  effluent 
fuel,  respectively.  The  treated  fuel  volumes  were  1.36  and  3.3  ml/g-ads  at  10  and  50  ppmw  S  in 
effluent  fuel,  respectively.  The  sample  Ads-2  has  a  good  level  of  surface  area,  total  metal  oxide, 
as  well  as  large  pore  size.  The  sample  Ads-1  has  lower  total  weight  percentage  of  metal  oxide, 
although  its  surface  area  is  slightly  higher  than  that  of  Ads-2.  The  sample  Ads-3  has  both  lower 
total  weight  percentage  of  metal  oxide  and  surface  area  than  that  of  Ads-2.  The  sample  Ads-4 
has  the  lowest  surface  area  even  if  its  metal  oxide  weight  percentage  is  high.  The  sample  Ads-5 
has  the  smallest  pore  size  even  if  there  is  a  sufficiently  high  surface  area.  Sample  6  has  similar 
weight  percentage  of  metal  oxide  and  a  slightly  lower  surface  area,  and  its  support  material  has 
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different  components  ratio;  therefore,  its  desulfurization  performance  is  also  lower  than  that  of 
sample  Ads-2.  Obviously  the  significant  low  surface  area  contributes  to  the  low  desulfurization 
capacity  of  sample  Ads-4.  As  a  conclusion,  with  metal  oxide  weight  percentage,  surface  area, 
and  pore  size  being  all  high,  the  desulfurization  performance  could  reach  the  best,  as  see  for 
sample  Ads-2.  The  XRD  of  Ads-2  is  shown  in  Figure  22.  The  XRD  graph  was  amorphous  with 
broad  peaks. 


Data*  01/29/14  08i44  S tap  t  0.020°  Cnt  Tima:  0.600  Sao . 
Rang*:  10.00  -  70.00  (Dag)  Cont.  Scan  Rata  :  2.00  Dag/min. 


Dag. 

Figure  22  XRD  results  ofNiQ-Ce02/Al20s-Si02  adsorbent 


3.4.5  Adsorbent  regeneration  and  capacity  recovery 

Regeneration  of  used  adsorbents  is  an  important  practical  issue,  since  the  sulfur  adsorption 
capacity  is  not  very  high  due  to  the  fact  that  the  sulfur  concentration  in  jet  fuel  is  high  and  also 
the  required  sulfur  concentration  in  treated  fuel  is  very  low.  For  regeneration  study  to  the 
currently  optimized  adsorbent,  sulfur  saturated  adsorbent  (used  previously)  was  heated  at  250  °C 
for  30  min  in  flowing  helium  for  removing/evaporating  fuel  held  in  pores,  and  then  heated  at  a 
raised  temperature  for  2  h  in  helium  to  regenerate.  The  four  different  raised  regeneration 
temperatures  were  300,  400,  500  and  600  °C.  The  desulfurization  performances  of  regenerated 
adsorbents  due  to  different  regeneration  temperatures  are  shown  in  Fig.  23.  Adsorbent 
regenerated  at  500  and  600  °C  had  better  desulfurization  performance  than  those  regenerated  at 
300  and  400  °C,  and  at  low  sulfur  concentration  level  in  effluent  fuel  the  adsorbent  regenerated 
at  500  °C  had  the  best  performance  with  the  high  volume  of  fuel  treatment  per  gram  of  adsorbent. 
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Treated  fuel  (ml/g~ads) 


Figure  23  Adsorbent  regeneration  at  different  temperatures 
fCo:  1140  ppmw;  LHSV:  0.63  h'1 ;  particle  size:  0-125  pm) 

After  the  first  regeneration  test,  multiple  cycles  of  regeneration  tests  were  performed  by  the 
repeated  treatment  of  the  saturated  adsorbent  at  500  °C  after  each  test  in  a  cycle.  The  four-cycle 
breakthrough  curves  and  the  adsorption  capacities  corresponding  to  sulfur  content  in  effluent  fuel 
are  shown  in  Fig.  24  and  Fig.  25.  The  sulfur  adsorption  capacity  of  the  regenerated  adsorbent 
decreased  with  the  increase  of  regeneration  cycles.  At  a  breakpoint  of  50  ppmw'  S  in  effluent  fuel, 
the  sulfur  adsorption  capacity  was  2.95  mg-S/g-ads  for  fresh  adsorbent,  and  was  2.21,  1.95,  1.76 
and  1.49  mg-S/g-ads  for  cycles  1  through  4.  In  cycle  1  through  4,  about  75%,  66%,  59%  and 
51%  sulfur  adsorption  capability  wras  recovered  in  the  regenerated  adsorbent.  The  measured  BET 
surface  areas  are  256,  253,  250,  222  and  219  m2/g  respectively  for  the  fresh  and  regenerated 
adsorbents  through  cycles  1  to  4.  The  BJH  average  pore  diameter  was  around  91±1  A  for  all 
adsorbents.  The  pore  size  distributions  are  shown  in  Fig.  26.  The  reduction  of  adsorption 
capacity  is  strongly  correlated  to  the  decrease  of  the  surface  area,  which  reflects  the  decrease  of 
the  active  sites  in  the  adsorbent  for  sulfur  adsorption  [54]. 


Figure  24  Four-cycle  breakthrough  curves 
(Co:  1140 ppmw;  LHSV:  0.63  h'1 ;  particle  size:  0-125  pm) 
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Sulfur  content  in  effluent  fuel  (ppmw) 

Figure  25  Capacity  versus  sulfur  content  in  effluent  fuel  for  four-cycle  tests 
(Cq:  1140 ppmw;  LHSV:  0.63  Kl;  particle  size:  0-125  pm) 


Figure  26  Pore  size  distribution  of  regenerated  adsorbent 
(Unit  of  pore  size  W:  A) 

4.0  Autothermal  reforming  of  n-dodecane  (surrogate  of  JP-8)  and  Jet-A  fuels 

The  conversion  of  hydrocarbon  fuel  to  hydrogen  or  syngas  can  be  made  based  on  three  main 
mechanisms.  These  are  steam  reforming  (SR),  partial  oxidation  (POX),  and  autothermal 
reforming  (ATR)  [55].  These  technologies  are  distinguished  by  whether  steam  or  oxygen,  or  a 
mixture  of  steam  and  oxygen,  is  used.  Production  of  syngas  in  the  current  industry  is  dominated 
by  a  steam  reforming  method  [56,57].  The  three  fuel  reforming  reactions  are  expressed  through 
Eqs.  (19)-(20). 

Steam  reforming:  CmH„  +  m\l20=mC0+(m+\n)\l2,  AH  >  0,  endothermic  (19) 

Partial  oxidation:  CmIln  +-j/w02=/wC0+(-f  n)H,,  AH  <  0,  exothermic  (20) 

Autothermal  reforming: 

CmHn  +  \mW20  +  \m02=mC0+{)2m+\n)Y\2,  A//<0,  exothermic  (21) 
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Steam  reforming  (SR)  is  a  high  endothermic  reaction  and  requires  a  lot  of  heat  from  an 
external  source.  It  is  an  approach  with  high  hydrogen  concentrations  and  higher  system 
efficiencies,  except  the  startup  time  is  typically  long  [58].  The  processing  system  tends  to  be 
heavy  for  steam  reforming  and  is  more  suitable  for  continuous  operation  under  a  steady  state, 
rather  than  for  onboard  fuel  cell  stacks  with  frequent  load  variations  [59],  In  a  partial  oxidation 
(POX)  process,  hydrocarbon  fuels  are  converted  to  syngas  by  partially  oxidizing  with  oxygen. 
This  is  a  high  exothermic  reaction  process,  and  generally  the  operating  temperature  ranges  from 
1100  °C  to  1200  °C  in  order  to  prevent  coking  in  the  reactor.  Compared  with  steam  reforming, 
partial  oxidation  has  the  advantages  of  compactness,  fast  start-up,  and  rapid  responses 
[58,60,61].  Depending  on  the  presence  or  absence  of  a  catalyst,  POX  has  two  types: 
homogeneous  POX  or  heterogeneous  catalytic  POX.  Homogeneous  POX  is  the  reactions  of  fuels 
with  the  oxygen  in  air  at  high  temperature  and  high  pressure  in  the  absence  of  any  catalyst  for 
producing  syngas  [1].  The  main  advantage  is  its  compatibility  with  various  fuels  from  natural  gas 
to  gasified  coal,  which  may  have  low  hazard  emissions  such  as  NOx  or  SOx.  However,  the 
hydrogen  production  efficiency  of  homogeneous  POX  is  relatively  low  due  to  the  very  high 
temperature  required  to  partially  oxidize  fuels  [62].  For  instance,  the  noncatalytic  process  for 
gasoline  reforming  requires  temperatures  in  excess  of  1000  °C  [63].  The  presence  of  a  suitable 
catalyst  can  reduce  the  operating  temperature  to  800-900  °C,  which  enables  the  application  of 
common  materials  such  as  stainless  steel  to  construct  a  reactor  and  also  increases  system 
efficiency  [64],  Regularly,  the  catalyst  can  only  be  employed  if  the  sulfur  content  in  the  fuel  is 
below  50  ppm  in  order  to  avoid  catalyst  poisoning  [65,66],  But  some  recent  studies  also  show 
POX  catalysts  can  tolerate  sulfur  of  a  content  of  up  to  a  couple  of  hundreds  ppm  [67-69], 
Autothermal  reforming  (ATR)  is  a  combination  of  exothermic  partial  oxidation  sequentially 
followed  by  endothermic  steam  reforming.  The  heat  released  by  POX  can  sustain  SR  and  the 
overall  ATR  reaction  is  thermally  neutral  or  slightly  exothermic.  The  operating  temperature  is 
usually  in  the  range  of  900  °C  to  1150  °C  in  reactors  with  a  lower  pressure  compared  to  POX.  So 
far,  reformers  with  pressure  between  1  and  80  bar  have  been  designed  and  built  [1].  The  H2  to 
CO  molar  ratio  in  ATR  reformate  is  about  2,  which  is  more  favorable  than  the  ratio  in  POX.  An 
external  heat  source  is  not  required  for  a  steady-state  ATR  process.  However,  startup  of  the 
oxidation  does  not  occur  at  ambient  conditions  and  requires  some  energy.  The  startup  phase  of 
ATR  is  called  light-off,  which  is  correspond  to  about  10%  total  oxidation  conversion  of  the  fuel 
in  a  typical  case  and  is  characterized  by  the  light-off  temperature.  The  light-off  temperature  is 
affected  by  fuel  type,  O2/C  ratio,  and  types  of  catalyst.  For  all  fuels,  noble  metal  has  higher 
activities  than  Ni  and  leads  to  relatively  lower  light-off  temperatures  [70],  The  light  off 
temperature  of  diesel  was  reported  by  Kang  et  al.  [71]  as  250  °C  and  by  Borup  et  al.  [72]  as  270 
°C.  The  light  off  temperature  of  Jet  A-l  fuel  was  reported  by  Karatzas  et  al.  [19]  as  300  °C.  The 
advantages  and  disadvantages  of  these  three  reforming  technologies  are  summarized  and 
compared  in  Table. 


Table  5  Comparison  of  reforming  technologies^ 


Technologies 

Advantages 

Disadvantages 

Steam  reforming 

Most  industrial  experience 

Highest  air  emissions 

Oxygen  not  required 

Heavy  system 

Lowest  temperature 

Highest  H2/CO  ratio 

Slow  startup 

Heat  source  required 

28  |  P  a  g  e 


Project  final  report  August  5th,  2015 


Partial  oxidation 


Higher  sulfur  tolerance 
No  heat  source  required 
Compact  system 
Fast  startup 


Low  H2/CO  ratio 
Highest  temperature 
coke  formation 
Oxygen  or  air  required 
Too  much  heat  produced 


Autothermal  reforming  Medium  temperature  Limited  experience 

No  heat  source  required  Oxygen  or  air  required 
Favorable  H2/CO  ratio 
Least  coke  formation 
Relatively  compactness 


Considering  the  requirements  of  size  and  weight  of  a  fuel  processing  system,  rapid  startup  and 
dynamic  responses  [73],  and  the  converting  ability'  of  complex  fuels  such  as  diesel  and  jet  fuels, 
the  ATR  approach  has  been  chosen  to  be  the  best  solution  for  onboard  reforming  applications 
[74,75].  Also,  because  high-temperature  SOFC  is  selected  to  be  the  end-user  unit,  which  can 
directly  consume  carbon  monoxide  the  same  as  hydrogen,  no  post-processing  reactors  such  as  a 
water  gas  shift  (WGS)  reactor  and  preferential  oxidation  (Prox)  reactor  is  required.  This  reduces 
the  system  complexity  and  the  difficulty  of  maintenance.  The  purpose  of  both  WGS  and  Prox  is 
to  selectively  decrease  carbon  monoxide  concentration  in  the  syngas  reformate  to  a  very  low 
level  of  1%  or  below  for  a  high  temperature  PEFC,  and  less  than  50  ppm  for  a  regular  PEFC 
[76,77]. 

Ni-based  catalysts  for  reforming  have  been  widely  used  for  many  years  because  of  their 
activity  and  low  cost.  The  first  catalyst  with  Ni  as  the  active  metal  for  reforming  was  reported  by 
Prettre  and  his  coworkers  [78].  Then  a  lot  of  research  works  on  Ni  based  catalyst  for  reforming 
were  carried  out.  However.  Ni  catalysts  have  inherent  challenges  such  as  sulfur  poisoning, 
carbon  formation  and  sintering  [20].  Thus  noble  metals  were  introduced  in  catalysts  for 
reforming  reactions  [79,80]  but  the  high  price  of  noble  metals  is  a  challenge.  Recently,  the  study 
of  bimetallic  metal  catalysts  by  introducing  noble  metals  into  Ni-based  catalysts  has  become  a 
popular  approach.  It  is  believed  that  the  bimetallic  metal  catalyst  attains  both  benefits  of  nickel 
and  noble  metal  [81],  and  improves  catalyst  performance  in  sintering  resistance  and  even 
distribution  of  temperatures  [82,83].  In  noble  metals,  Pt  has  a  high  activity  for  oxidation  but  low- 
activity  for  steam  reforming.  Pd  has  a  better  steam  reforming  activity  than  Pt  but  is  sensitive  to 
coke  formation.  Rh  and  Ru  are  very  good  catalyst  candidates  because  of  their  good  activity  for 
both  oxidation  and  steam  reforming  reactions.  The  price  of  Rh  is  relatively  high  for  noble  metals, 
and  Ru  is  relatively  cheaper. 

Autothermal  reforming  of  jet  fuels  has  been  investigated  by  several  research  teams.  Lenz  et 
al.  [84]  studied  autothermal  reforming  of  desulfurized  Jet  A-l  with  a  15  kWt  test  rig  and  the  best 
energy  conversion  efficiency  of  78.5%  was  obtained  at  air  to  fuel  ratio  of  0.28  and  steam  to 
carbon  ratio  of  1.5.  Pasel  et  al.  [85,  86]  studied  autothermal  reforming  of  commercial  Jet  A-l  on 
a  5  kWe  scale.  For  a  C13-C15  alkane  mixture  surrogate  fuel  they  reported  an  optimized  energy- 
conversion  efficiency  of  80%  at  oxygen  to  carbon  molar  ratio  O2/C  =  0.43  and  steam  to  carbon 
molar  ratio  S/C  =  1.9.  Karatzas  et  al.  [19]  tested  autothermal  reforming  of  commercial  Jet  A-l 
and  the  optimized  energy  conversion  efficiency  was  only  42%,  which  led  to  the  conclusion  that 
the  high  sulfur  content  in  the  fuel  had  a  detrimental  effect  on  the  reformer  performance. 
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There  are  several  issues  to  be  addressed  for  the  design  and  operation  of  a  reforming  system 
for  heavy  hydrocarbon  fuels  [87],  First,  local  hot-spots  may  exist  in  the  reforming  catalyst  due  to 
the  non-uniform  temperature  distribution,  and  local  catalyst  deactivation  can  be  caused  by  the 
local  high  temperature.  Second,  the  large  carbon  contents  may  cause  coking,  which  can 
significantly  decrease  the  effectiveness  of  the  catalyst  and  thus  reforming  efficiency  [88]  Hard 
coke  such  as  graphitic  is  unreactive  with  hydrogen  and  can  block  active  sites  of  the  catalyst. 
Third,  inhomogeneous  mixing  of  reactants,  air,  fuel  and  steam,  can  cause  both  local  hot-spots 
and  coke  formation  in  the  catalyst  surface.  Fourth,  as  discussed  in  the  previous  section,  sulfur 
components  could  deactivate  the  reforming  catalyst  and  poison  electrodes  in  fuel  cells.  The 
catalyst  deactivation  can  be  caused  by  the  formation  of  strong  metal-S  bonds  which  modify  the 
surface  chemistry.  The  chemisorbed  sulfur  onto  the  active  catalyst  sites  also  prevents  reactant 
access  [89,90].  Generally,  coke  formation  can  be  suppressed  by  controlling  the  fuel  evaporation 
temperature  [91]  and  optimizing  the  mixing  ratio  of  fuel,  water  steam,  and  air.  The  catalyst 
poisoning  problem  due  to  sulfur  can  be  solved  by  fuel  desulfurization  technology  [92-94], 

The  autothermal  reforming  reactor  consists  of  a  thermal  zone  where  partial  oxidation 
happens,  which  generates  heat  to  drive  the  steam  reforming  reactions  in  a  downstream  catalytic 
zone  [95,96].  The  thermal  zone  may  take  20%  of  the  top  catalyst  bed  for  catalytic  POX. 
Therefore,  the  temperature  distribution  in  the  axial  direction  of  the  reformer  is  non-uniform.  The 
temperature  profile  has  a  sharp  rise  to  the  peak  in  the  POX  zone  and  then  a  decrease  due  to  the 
endothermic  reactions  to  a  relatively  low  and  flat  level  in  the  SR  zone  [97,98].  The  non-uniform 
axial  temperature  distribution  could  cause  the  problem  of  so  called  “hot-spots”,  which  induce  the 
potential  risk  of  local  catalyst  deactivation  due  to  thermal-induced  mechanisms  such  as  sintering. 
Sintering  is  defined  as  a  thermal  process  that  produces  a  decline  in  surface  area  of  the  active 
ingredient  or  the  support  of  the  catalyst,  which  results  in  a  decline  in  the  observed  rate  of 
reactions  [99].  As  Qi  [21]  reported  in  the  ATR  of  gasoline,  due  to  the  temperature  gradient 
greater  than  150  °C  in  the  axial  direction,  the  catalyst  rapidly  deactivated  and  the  gasoline 
conversion  decreased  from  100%  to  95%  after  40  hr  operation. 

The  primary  measure  to  minimize  hot-spots  is  to  employ  proper  reactor  materials,  catalyst 
support  structures,  and  flow  configurations,  which  favor  effective  heat  transfer  and  more  uniform 
axial  temperature  distribution.  Stainless  steel  can  be  used  to  construct  the  reactor  because  of  its 
relatively  good  thermal  conductivity  and  high  temperature  tolerance.  A  traditional  pellet  catalyst 
with  a  packed-bed  configuration  has  poor  heat  transfer  performance  [100]  and  it  has  been 
suggested  that  it  be  replaced  by  metallic  monoliths,  foams,  wire-gauzes,  or  microchannel 
reactors  [101,102],  Flow  with  high  turbulence  can  improve  hot-spots  since  the  turbulent  flow 
enhances  the  heat  transfer  coefficient  between  the  flow  and  solid  walls.  Therefore,  structured 
flow  paths  allowing  high  flow  rates  are  preferred.  The  selection  of  appropriate  channel  diameters 
and  geometries  of  the  catalyst  monolith  is  to  achieve  effective  mixing,  and  high  Sherwood  and 
Nusselt  numbers,  which  are  equivalent  to  high  mass  and  heat  transfer  coefficients.  Other  issues, 
such  as  pressure  drops,  also  need  to  be  considered  to  be  within  acceptable  limits. 

Coke  is  a  high-molecular-weight  polymer  with  low  hydrogen  content.  Formation  of  coke  on 
the  catalyst  surface  is  thought  to  occur  in  many  instances  by  polymerization  of  aromatic 
compounds  originally  present  or  formed  in  reactions  [99],  Olefin  and  aromatic  contents  in  diesel 
and  jet  fuels  are  precursors  of  coke  formation  [98,103],  Coke  formation  could  cause  degradation 
of  the  reformer  performance  and  reduces  its  life  time  significantly  [104-106],  Souza  et  al.  [107] 
reported  that  in  the  ATR  of  methane  over  Pt/Zr02  catalyst,  the  fuel  conversion  decreased  from 
80%  to  65%  after  30  h  of  operation  in  the  existence  of  coke  formation.  Yoon  et  al.  [108]  reported 
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that  in  the  ATR  of  synthetic  diesel,  due  to  coke  formation  the  fuel  conversion  decreased  from 
100%  to  90%  and  the  reformer  efficiency  decreased  from  65%  to  45%  after  40  h  of  operation. 
High  temperature  is  required  to  minimize  coke  formation.  Holladay  [109]  summarized  the 
minimum  reaction  temperatures  required  for  avoiding  coke  formation  during  isooctane  reforming 
at  thermodynamic  equilibrium.  For  C>2/C=0.5  and  H20/C=l,  the  minimum  required  temperature 
is  1030  °C  and  for  OfC=\  and  H2O/O2,  the  temperature  is  drastically  reduced  to  575  °C.  The 
required  temperature  for  ATR  is  higher  than  that  for  SR  [97,109,110],  but  it  was  proved  that 
ATR  produces  less  coke  compared  with  SR  and  POX  [111,112].  Promoted  catalysts  have  a 
positive  effect  in  coking  resistance,  and  generally,  Rh-based  catalysts  give  better  performance 
than  Ni-based  catalysts  in  coking  prevention. 

Ethylene  (C2H4)  is  reported  to  be  the  main  reason  for  rapid  carbon  formation  in  the  reforming 
process  [98]  A  large  quantity  of  ethylene  can  be  produced  by  pyrolysis  of  heavier  hydrocarbons 
at  the  local  fuel-rich  zones  if  the  mixing  of  fuel,  air,  and  steam  is  not  homogeneous.  Ethylene  can 
be  decomposed  into  carbon  in  the  absence  of  oxygen  and  water.  In  order  to  investigate  the 
mechanism  of  carbon  formation  caused  by  ethylene  pyrolysis,  Yoon  et  al.  [108]  designed  a  diesel 
ATR  processing  system  to  study  the  correlation  between  carbon  formation  and  O2/C  and  H2O/C 
ratios  under  the  reforming  conditions  of  SR,  POX,  and  ATR.  A  blank  reactor  and  a  catalyst- 
loaded  reactor  were  used  to  distinguish  the  carbon  formation  performance  between  homogeneous 
(in  a  blank  reactor)  and  heterogeneous  (in  a  catalyst-loaded  reactor)  reactions.  Their  testing 
results  for  two  reactors  under  the  ATR  condition  show  that  the  homogeneous  reaction  produced 
much  more  ethylene  than  the  heterogeneous  reaction,  and  therefore  they  concluded  that  almost 
all  the  ethylene  is  produced  in  a  homogeneous  reaction  at  the  reactor  entrance  for  ATR  reactors. 
The  comparison  of  CO  level  in  the  reformate  shows  that  both  POX  and  SR  have  more  carbon 
deposition  than  ATR.  Figure  27  shows  the  detected  carbon  by  TPO  (temperature  programmed 
oxidation  method)  analysis  for  three  reforming  approaches,  and  the  two  SR  lines  represent  SR 
with  different  II2O/C  ratios.  Apparently,  insufficient  steam  resulted  in  reforming  performance 
degradation,  and  carbon  formation  is  prone  to  happen  in  the  absence  of  oxygen.  In  the  real  fuel 
ATR  process,  most  oxygen  is  consumed  by  fuel  decomposition  in  the  homogeneous  reaction  at 
the  entrance.  A  lack  of  oxygen  in  the  downstream  section  causes  formation  of  ethylene.  Because 
ethylene  is  practically  exposed  to  SR  in  the  catalyst  zone,  it  causes  carbon  deposition.  As  a 
conclusion,  the  control  of  H2O/C  and  O2/C  ratios  of  the  feeding  gases  plays  an  essential  role  in 
suppression  of  carbon  formation. 


Figure  27  Carbon  detected  by  TPO  for  C2H4  reforming 
(TPO:  air  500ml/min,  18-900  °C,  10  °C  /min)  [108] 
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Complete  fuel  evaporation  and  homogeneous  mixing  with  air  and  steam  are  a  big  challenge  in 
the  ATR  reactor  design  for  long-chain  hydrocarbon  fuels.  Many  undesired  problems  can  result  in 
insufficient  fuel  evaporation  or  inhomogeneous  mixing  [113].  First,  coke  formation  can  be 
formed  immediately  if  liquid  fuel  contacts  the  catalyst  due  to  the  presence  of  aromatics  in  diesel 
and  jet  fuels.  Second,  unexpected  hot-spots  can  be  caused  by  the  local  occurrence  of  insufficient 
steam  or  excess  oxygen  if  the  mixing  of  fuel,  oxygen,  and  steam  is  not  homogeneous  [114]. 
Stable  and  sustainable  hydrogen  throughput  also  requires  homogeneous  and  constant  mixing. 
Lindstrom  [115]  compared  the  performances  of  three  reactors  with  different  mixing  chamber 
designs.  In  the  M2  reactor  the  diesel-slip  (diesel  not  converted)  in  the  reformate  was  1,500  ppm, 
and  the  H2  volumetric  percentage  in  the  production  suddenly  decreased  from  35%  to  25%  after 
200  min  of  operation.  In  the  M4  reactor  the  diesel-slip  was  less  than  15  ppm,  and  the  H2  was 
stable  at  35%  level  after  450  min  of  operation.  In  the  M5  reactor  with  an  improved  fuel  injection 
system,  the  diesel-slip  was  controlled  less  than  10  ppm  and  the  H2  was  stable  at  40%  level. 
Another  issue  related  to  complete  fuel  evaporation  and  homogeneous  mixing  is  a  safety  concern. 
Since  the  boiling  temperature  for  diesel  and  jet  fuels  is  higher  than  the  auto-ignition  temperature, 
if  the  evaporation  and  mixing  do  not  complete  rapidly  and  uniformly,  the  formation  of  local 
oxygen-rich  zones  is  possible.  In  these  zones,  the  volume  ratio  of  fuel  and  oxygen  can  be  lower 
than  the  critical  values  and  may  end  up  in  an  explosion,  which  is  quite  dangerous  in  a  well-sealed 
space. 


4.1  Theory  and  thermodynamics  analysis  of  autothermal  reforming 

The  following  thermodynamics  analysis  for  a  reforming  process  will  make  clear  of  the 
required  proper  oxygen  to  carbon  ratio  and  the  range  of  reaction  temperature  that  an  autothermal 
reaction  can  sustain  without  external  heating  and  cooling.  For  a  general  hydrocarbon  or 
oxygenate  fuel,  CmHn07„  the  autothermal  reforming  reaction  stoichiometry  can  be  expressed  as 

CmH„0.  +  a02  +  M  I20=cCH4  +  dCO  +  eC02+  fii2  (22) 

where  CO2  comes  from  the  water-gas-shift  reaction  and  CH4  from  methanation  reaction  as 
follows: 

C0+H20=C02+H2  (23) 

CO+3H2=CH4  +H20  (24) 

C02  +  4H2=CH4  +  2H20  (25) 

Based  on  atomic  balance,  the  stoichiometric  coefficients  can  be  expressed  as  follows: 

d  =  2m  —  z  — 2a  — b  — 2c  (26) 

e ——m  + z  +  2a  +  b  +  c  (27) 

f  =  ^  +  b-2c  (28) 

Since  both  hydrogen  and  carbon  monoxide  are  fuels  for  SOFC,  the  total  yield  is 

d+ f  =  2m  +  ^-z-2a-4c  (29) 

which  implies  that  the  theoretical  maximum  yield  can  be  achieved  when  no  methane  is  produced. 
For  each  of  the  WGS  and  methanation  reactions,  at  equilibrium  state  the  system  Gibbs  energy  is 
zero  and  the  equilibrium  equation  is  satisfied, 

-AG°  =  RT\nK  (30) 
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The  relation  of  equilibrium  constant  K  and  temperature  can  be  plotted  as  shown  in  Fig.  28. 
Commonly  the  working  temperature  of  the  ATR  processor  is  in  the  range  of  650  °C  to  900  °C 
[1 16].  In  this  temperature  range  the  equilibrium  constant  K  for  reactions  by  Eq.  (24)  and  (25)  are 
almost  zero,  which  indicates  that  production  of  methane  due  to  methanation  reactions  is 
ignorable.  Nevertheless,  both  CO  and  CO2  exist  in  the  products  of  autothermal  reforming. 


Figure  28  Relation  of  equilibrium  constant  K  and  temperature 


Rase  [99]  proposed  a  simple  method  to  estimate  the  equilibrium  constant  value  for  ideal  gases 
in  a  reaction.  Since  the  standard  free  energy  can  be  expressed  by  standard  enthalpy  and  entropy, 

AG°  =AH-TAS°  (31) 

Dividing  both  sides  by  T,  the  differential  term  can  be  obtained  as 

d(AG"  IT)  \(dAH\  AH  dAS° 


dT 


dT  , 


dT 


(32) 


Since  dHj=TdSj+VdP ,  differentiating  the  equation  with  dT  and  putting  the  result  back  into  Eq. 
(32),  the  following  equation  is  obtained: 

d(AG°!T) 


Combining  Eqs.  (30)  and  (33)  gives 


dT  T2 

d\nK  AH 


(33) 


(34) 


dT  RT2 

The  AH  can  be  obtained  through  integrating  the  standard  empirical  heat  capacity  as  shown  in  Eq. 
(35). 


,(T) 


(35) 


where  P  represents  for  products  and  R  represents  for  reactants  and  v  is  the  stoichiometric 
coefficient.  CpfT)  is  the  heat  capacity  of  species  j  as  a  function  of  temperature.  £engel  [117]  and 
Chase  [118]  summarized  the  empirical  expression  of  CpfT)  as  follows, 

CPJ  (T)  =  Aj+  BT  +  CjT2  +  DT}  (36) 
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A  through  D  are  empirical  coefficients  for  heat  capacity.  For  the  substances  the  empirical 
expression  of  which  are  not  given,  Cpj(T)  can  be  calculated  by  Method  of  Joback  [119]  based  on 
group  contributions  analysis.  By  integrating  Eq.  (36),  the  system  enthalpy  change  is  given  as, 


ah=i„+±vjhj-±vjhj 

(37) 

P  R 

BT 2  CT3  Z)T4 

H  =  AT +  +  +  ' 

'  7  2  3  4 

(38) 

The  right-hand-side  term  of  Eq.  (34)  can  be  substituted,  and  after  integration,  Eq.  (34)  is 
expressed  in  a  new  form  as  given  by  Eq.  (39): 


In  AT  =  ■ 


R 


-kr+l'+h'fo-isj  K, 

1  P  R 


B,T  C,T 2 

K,=Al\nT  +  —  +  ^—  + 

J  J  2  6 


DjT 3 
12 


(39) 

(40) 


where  7//  and  Ik  are  integration  constants  that  can  be  evaluated  based  on  the  substance 
information  at  298K.  Using  this  equation,  the  equilibrium  constant  AT  can  be  estimated  for 
different  reaction  temperatures  if  the  reaction  stoichiometries  are  determined.  Take  the  following 
reforming  reaction  of  dodecane,  for  example, 


C12H2t)  +  x02  +  (24  -  2x)  H20=12C02  +  (37  -  2x)  H2  (41) 


The  effects  of  temperature  on  the  s>  stem  enthalpy  change,  system  Gibbs  free  energy  change  and 
equilibrium  constant  AT  at  various  O2  to  carbon  ratios  at  the  operation  temperature  between  800 
and  1500  K  are  shown  in  Fig.  29  to  Fig.  31.  Figure  29  indicates  that  the  O2/C  ratio  has  to  be 
greater  than  0.3  to  meet  the  exothermic  reaction  requirement.  If  O2/C  ratio  is  less  than  0.3  the 
reaction  becomes  endothermic,  which  cannot  self-sustain.  Figure  30  shows  that  the  reaction  can 
go  forward  when  the  O2/C  ratio  is  above  0.20.  Figure  31  shows  the  operation  temperature 
doesn’t  influence  the  equilibrium  constant  much.  In  contrast,  high  oxygen  to  carbon  ratio  can 
increase  the  value  of  the  equilibrium  constant. 


Figure  29  System  enthalpy  change  versus  temperature 
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Figure  30  System  Gibbs  free  energy'  change  versus  temperature 


Figure  31  Equilibrium  constant  lnKp  versus  temperature 

4.2  Monolithic  reformer 

Monolithic  reactor  is  one  popular  type  of  reactor  in  reforming  systems.  Monoliths  are 
continuous  structures  with  well-defined  geometries  involving  parallel  and  identical  channels  with 
small  diameters.  The  shape  of  the  cross-section  of  the  channels  can  be  square,  sinusoidal, 
triangular,  hex,  round,  and  so  on  [120].  Fig.  32  illustrates  three  monoliths  with  different  cell 
geometries  [121].  Compared  with  a  packed  bed,  a  monolithic  or  honeycomb  structure  has  the 
advantages  of  low  pressure  drop  and  high  surface  area/volume  ratio.  The  large  open  front  area 
(OFA)  and  straight  parallel  channels  in  monoliths  ensure  small  flow  resistance,  even  at  high  flow 
rates.  The  large  surface  area  is  usually  achieved  by  depositing  a  high-surface-area  carrier,  a 
catalyst  that  can  be  impregnated  into  the  channel  walls.  The  deposited  catalyst  carrier  is  called 
washcoat,  and  its  thickness  is  determined  by  the  geometry  of  the  channel  and  coating  method. 
Too  thick  a  washcoat  could  result  in  an  increase  of  pressure  drop  to  an  unacceptable  level. 
Besides,  the  increase  in  cell  density  could  cause  a  significant  increase  in  both  surface  area  and 
pressure  drop  [120]. 
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Figure  32  Monoliths  with  different  cell  geometries  [121] 

Synthetic  cordierite  is  the  most  commonly  used  ceramic  for  monolithic  catalyst  support.  Its 
characteristics  of  low  thermal  expansion  coefficient  in  the  magnitude  of  10'6  K'1,  high 
mechanical  strength  of  3,000  pounds  per  square  inch,  and  high  melting  point  over  1,300  °C  make 
it  the  preferred  material  to  apply  to  an  exothermic  reaction.  However,  the  thermal  shock 
resistance  of  the  monolith  can  be  influenced  by  the  washcoat,  which  usually  has  a  larger  thermal 
expansion  coefficient.  More  attention  needs  to  be  paid  to  this  thermal  shock  resistance  difference 
in  rapid  temperature  change  situations  such  as  startup  and  shutdown.  Particle  size  of  the  carrier 
and  thickness  of  the  washcoat  are  two  key  parameters  that  can  be  optimized  to  decrease  an 
undesired  impact.  One  drawback  of  a  cordierite  monolith  is  it  is  less  suitable  for  isothermal 
operations  due  to  its  low  thermal  conductivity.  High  temperature  resistant  aluminum-containing- 
steel  is  another  popular  material  for  monoliths.  Its  main  advantage  compared  to  cordierite  is  its 
potential  to  be  manufactured  with  thinner  walls,  which  leads  to  higher  cell  densities  and  lower 
pressure  drop.  The  15-20  times  higher  thermal  conductivity  than  ceramic  is  another  advantage 
that  makes  the  isothermal  operation  possible,  which  favors  fast  startup.  However,  its  application 
is  limited  to  electrically  heated  catalytic  converters  as  it  is  electrically  conductive  [120], 

The  most  important  characteristics  to  justify-  the  performance  of  monoliths  are  cell  density  ( n ), 
geometric  surface  area  (GSA),  open  frontal  area,  hydraulic  diameter  {D[),  bulk  density,  the 
thermal  integrity  factor  (TIF),  the  mechanical  integrity  factor  (MIF),  resistance  to  flow  (RJ),  bulk 
heat  transfer  ( Hs ),  and  light-off  (LOF)  [120]  The  asymptotic  Sherwood  ( Sh )  number  and  Nusselt 
( Nu )  number  are  two  dimensionless  numbers  that  are  used  to  characterize  heat  and  mass  transfer 
coefficients.  Both  numbers  have  constant  values  for  fully  developed  laminar  flow-  away  from  the 
entrance  in  channels  with  fixed  diameter  and  shape.  Also,  the  flow  in  monolithic  channels  is 
always  laminar  since  the  diameter  is  quite  small  and  the  Reynolds  number  is  low.  Pressure  drop 
across  the  substrate  depends  linearly  on  flow  velocity  and  length  but  inversely  on  the  square  of 
the  hydraulic  diameter.  Sherwood  number,  Nusselt  number,  and  pressure  drop  are  defined  as 
follows,  respectively: 


Sh  =  K-Dh/dAB 

(42) 

Nu  =  h-  Dhl  kf 

(43) 

A__  CfVJL 

P  ADl(  OFA) 

(44) 

where  K  is  the  mass  transfer  coefficient,  dAB  is  the  mass  diffusivity,  h  is  the  heat  transfer 
coefficient,  k y  is  the  thermal  conductivity  of  fluid,  C/  is  the  friction  coefficient,  A  is  the  cross 
section  area  of  the  substrate,  L  is  the  length  of  the  substrate,  and  Ve  is  volumetric  flow  rate. 
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For  cordierite  monoliths,  thin  and  ultrathin  wall  structures  with  high  cell  density  have  been 
investigated  to  minimize  thermal  mass  and  maximize  surface  area  [122,  123].  It  has  been  found 
that  the  thin  wall  substrates  have  40%  lower  heat  capacity,  50%  lower  mass,  and  60%  higher 
GSA  than  standard  substrates.  However,  higher  pressure  loss  is  also  observed  because  of  denser 
cells  and  increased  flow  resistance.  The  performances  of  square,  triangular,  and  hexagonal  (hex) 
cells  have  also  been  compared  to  study  the  effect  of  different  cell  shapes.  The  square  and 
triangular  cells  have  been  proven  to  be  the  most  cost  effective  in  terms  of  extrusion  die  cost 
[124].  Experimental  results  show7  that  a  hex  cell  has  7%  lower  thermal  mass  than  a  square  cell, 
and  the  triangular  cell  has  a  13%  higher  thermal  mass  than  a  square  cell.  However,  by  including 
the  effects  of  the  heat  transfer  factor  and  GSA  on  a  light-off  factor  (LOF)  and  conversion 
efficiency  factor  (CEF).  the  square  cell  offers  an  equivalent  or  even  slightly  better  performance 
than  that  of  hex  cells  because  of  the  inherent  high  GSA  of  a  square  cell,  although  the  thermal 
mass  for  a  square  cell  is  higher.  On  the  other  hand,  the  LOF  and  CEF  values  of  a  triangular  cell 
appear  to  be  even  higher  than  those  of  a  square  cell.  Considering  that  in  real  cases  there  is  little 
or  no  flow  in  the  acute  comer  regions,  the  effective  performance  of  a  triangular  cell  is  not  as 
good  as  a  square  cell.  Furthermore,  the  pressure  drop  of  a  triangular  cell  substrate  is  nearly  30% 
higher  than  that  of  a  square  cell  substrate,  which  makes  a  square  cell  structure  the  preferred 
choice  in  terms  of  overall  performance.  A  hex  cell  substrate  has  approximately  10-12%  less 
pressure  loss  than  a  square  cell  substrate,  which  is  the  only  advantage  of  a  hex  cell  over  a  square 
cell  [120]. 

The  challenges  of  monolith  reactors  are  mainly  related  to  even  flow  distribution  and 
replacement  of  the  catalyst  [113].  A  uniform  distribution  is  essential  for  keeping  a  narrow  RTD, 
which  is  important  in  the  ATR  process  to  make  sure  that  a  partial  oxidation  reaction  does  not 
transform  into  a  total  oxidation  reaction.  A  supplemental  device  for  distributing  reactant  flow  is 
required,  and  the  system  becomes  more  complex.  Catalyst  deactivation  is  also  a  serious  problem, 
which  a  monolith  reactor  cannot  avoid.  Theoretically,  the  replacement  can  be  done  by 
disintegration  of  the  washcoat  from  the  monolithic  support,  but  it  is  obviously  difficult  to 
execute  and  the  capital  cost  is  also  significant.  One  additional  problem  monolith  reactors  may 
encounter  is  washcoat  loss.  The  high  flow  velocities  and  rapid  temperature  changes  could  lead  to 
the  loss  of  bonding  between  the  washcoat  and  monolith  walls.  Research  teams  in  Argonne 
National  Laboratory  [73],  Royal  Military  College  of  Canada  [21,125,126],  Forschungszentrum 
Juelich  GmbH  in  Germany  [86,114,127],  and  KTH-RIT  in  Sweden  [19]  have  done  a  significant 
amount  of  research  of  ATR  in  monolithic  reactors.  For  diesel  ATR  in  monolithic  reactor, 
Shigarov  et  al.  [128]  tested  several  catalyst  composites  and  the  optimum  operating  conditions 
were  specified  as  O2/C  ratio  of  0. 5-0.6,  S/C  ratio  of  1.5-1. 7  and  inlet  mixture  temperature  of  300- 
400  °C.  At  these  conditions  the  yield  of  H2+CO  was  2.88  L/g  fuel,  the  hydrogen  yield  was  18 
mol/mol  fuel,  and  the  H2/CO  ratio  was  3.5  in  the  products.  Karatzas  [19]  studied  a  monolithic 
ATR  reformer  of  5  kWe  for  diesel  (~10  ppm  S)  and  jet  fuel  (-200  ppm  S)  reforming.  At 
optimized  conditions  of  O2/C  ratio  of  0.49,  S/C  ratio  of  2.5  and  fuel  inlet  temperature  of  60  °C, 
reformer  efficiency  of  81%  and  H2  selectivity  of  96%  were  established  for  diesel.  And  at  the 
same  operating  conditions,  the  reforming  efficiency  of  jet  fuel  was  only  42%.  Lenz  et  al.  [84] 
investigated  the  ATR  of  desulfurized  jet  fuel  in  a  monolithic  reactor.  At  S/C  ratio  of  1.5,  air  to 
fuel  ratio  of  0.28,  and  GHSV  (Gas  hourly  space  velocity)  of  50000  h’1,  the  best  efficiency  of 
78.5%  was  reached. 

Besides  the  traditional  honeycomb  monoliths  carriers,  foam  monoliths  have  been  considered 
as  an  alternative  catalyst  carrier  for  their  advantages  of  better  heat  and  mass  transfer  properties. 
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Their  mass  transfer  coefficients  are  between  those  of  packing  particles  and  honeycomb 
monoliths,  and  their  heat  transfer  coefficients  are  efficient  compared  to  particle  beds.  The  flow 
through  the  foam  structure  follows  the  same  convective  fluid  mechanics  as  in  packed  beds  and 
better  turbulence  is  obtained  compared  to  honeycomb  monoliths.  The  pressure  drop  in  a  foam 
structure  is  reduced  by  a  factor  of  10  compared  to  a  packed  bed  since  porosities  are  almost  twice 
as  large.  But  the  pressure  loss  is  higher  than  that  through  honeycomb  monoliths  [129,130], 


4.3  Catalyst  preparation  for  the  reformer 

The  NiO-Rh  bimetallic  catalysts  were  prepared  by  the  washcoating  method  [21,131,132], 
Support  substrates  were  cordierite  monoliths  with  400  cpsi  and  a  wall  thickness  of  0.25  mm, 
length  of  60  mm  and  diameter  of  40  mm.  Picture  and  XRD  graph  of  the  used  cordierite  monolith 
are  shown  in  Fig.  33.  Pseudo-boehmite  (70  wt.%  AI2O3)  was  coated  as  the  first  layer  in  the 
monoliths  to  increase  the  total  surface  area.  Cerium  was  added  to  increase  the  sulfur  resilience 
capability  of  the  catalyst,  as  explained  by  Argonne  National  Laboratory  that  Ce  could  form  a 
stable  sulfide  in  the  temperature  range  of  ATR  to  serve  as  a  sulfur  sink  [21].  Potassium  was 
proved  to  be  the  most  effective  promoter  in  suppressing  coke  formation  in  ATR.  The  addition  of 
K  not  only  lowers  the  coking  tendencies  of  alumina  supports,  but  also  prevents  nickel  from 
catalyzing  the  decomposition  reactions  [133],  Lanthanum  was  added  as  a  promoter  to  improve 
the  catalytic  activity  of  nickel.  In  reforming  reactions,  the  addition  of  a  small  amount  of  La  could 
possibly  increase  the  nickel  catalyst  activity  significantly  [134],  Nickel  oxide  and  rhodium  were 
the  major  catalytic  components  in  the  proposed  catalysts. 


(a) 


Figure  33  (a[  Picture  of  the  used_  cordierite  monolith  (b)  XRD  graph  of  the  cordierite  monolith 
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Regent  grade  rhodium  chloride  hydrate  RhCl3.3H20  (Rh  >  39.5  wt.%),  cerium  acetate  hydrate 
Ce(CH2C00H)3.5H20  and  lanthanum  acetate  hydrate  La(CH2COOH)3  5H2O  (from  Hangzhou 
Ocean  Chemical  Co.,  Ltd.,  China),  nickel  acetate  hydrate  Ni(CH2C00H)2.4H20  (from 
Sinopharm  Chemical  Reagent  Co.,  Ltd.,  China),  and  potassium  nitrate  KNO3  (from  Sigma- 
Aldrich)  were  used  as  sources  of  metals  in  the  catalysts.  Before  coating  the  cordierite  monoliths 
were  cleaned  in  10%  nitric  acid  and  then  washed  by  distilled  water.  Followed  by  drying  at  60  °C 
overnight  the  clean  monoliths  were  calcined  in  air  at  815  °C  for  2  h.  After  cooling  the  monoliths 
were  ready  for  washcoating.  The  first  step  was  to  impregnate  the  cordierite  monoliths  into 
prepared  Al-Ce  sol  solution  for  20  min.  The  surplus  solution  remaining  in  the  monolithic 
channels  was  evacuated  by  N2  flow.  After  dry  ing  at  60  °C  overnight  the  monoliths  were  then 
calcined  in  air  at  815  °C  for  2  h.  The  process  was  repeated  until  the  desired  amount  of  AI2O3- 
CeC>2  was  coated  in  the  monoliths.  The  same  procedure  was  then  applied  to  the  monoliths 
repeatedly  with  La-K  solution,  Ni  solution  and  Rh  solution. 
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Figure  34  Pore  size  distribution  of  monoliths  with  different  treatments 

(Unit  of  pore  size  W:  A) 

Micromeritics  TriStar  II  3020  was  used  to  measure  the  surface  area  and  pore  size  distribution 
of  the  cordierite  monoliths.  Figure  34  shows  the  pore  size  distribution  of  monolith  without  any 
treatment,  monolith  treated  with  acid  cleaning,  monolith  treated  with  calcination,  monolith 
treated  with  both  acid  cleaning  and  calcination,  and  monolith  coated  by  7  wt.%  AI2O3  and  2 
wt.%  Ce02.  The  corresponding  surface  areas  were  0.525,  6.273,  0.285,  1.417  and  5.332  m2/g 
respectively.  It  is  clear  that  the  acid  cleaning  can  remove  undesired  matters  from  the  monolith, 
and  the  alumina  coating  layer  can  significantly  increase  the  surface  area  of  the  monolith,  as  well 
as  increase  the  number  of  pores  in  the  size  range  between  10  and  100  nm.  Since  structured 
monolith  support  with  large  surface  area  and  small  pores  is  desired  in  the  catalyst  preparation 
process,  in  further  tests  all  cordierite  monolith  supports  were  firstly  washed  with  10%  nitric  acid 
followed  by  calcination  in  air  at  815  °C  for  2  h.  And  before  the  promoters  and  active  components 
were  added,  several  layers  of  alumnia  were  washcoated  onto  the  monolith. 

4.4  Autothermal  reforming  system  design 

Figure  35  indicates  the  lab  scale  experimental  set-up  of  the  2.5  kWt  ATR  system  [135].  The 
liquid  fuel  and  water  were  delivered  to  the  reactor  by  Flash  100  HPLC  (High-performance  liquid 
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chromatography)  pumps  at  constant  flow  rates.  Water  was  heated  to  superheated  steam  in  a 
Paragon  kiln  before  entering  the  reactor.  Liquid  fuel  and  compressed  air  were  firstly  pre-heated 
in  a  kiln  through  separate  tubing  paths  and  then  heated  in  heating  hoses  to  preset  temperatures. 
The  H900  series  heating  hoses  are  made  by  Hillesheim,  with  diameter  of  80  mm,  length  of  50  cm, 
and  maximum  temperature  of  450  °C.  The  heating  hoses  temperatures  were  controlled  by  HT  43 
PID  controllers. 

Furnace 


Figure  35  Experimental  set-up  of  the  lab-scale  ATR  system 


The  preheated  liquid  fuel  was  injected  into  the  reactor  through  a  stainless  steel  Steinen  mist- 
jet  misting  nozzle.  The  misting  nozzle  generated  hollow  cone  fuel  spray  with  a  spray  angle  of 
90°.  The  preheating  of  fuel  could  lower  the  fuel  viscosity  and  surface  tension  and  therefore 
benefit  the  fuel  transportation  and  injection  through  the  nozzle  [19].  Superheated  steam  and  air 
were  evenly  distributed  in  the  tubing  system  as  shown  in  Fig.  36  before  entering  the  ATR 
reformer.  Both  steam  and  air  were  injected  into  the  reactor  through  four  holes  in  the  wall.  The 
steam  injection  was  positioned  upstream  and  the  air  injection  was  positioned  downstream  to 
prevent  possible  fuel  ignition  [85]  All  the  tubing  was  wrapped  with  heat  insulation  materials  to 
reduce  heat  loss.  Followed  by  injection  the  fuel  spray  was  firstly  evaporated  by  the  heat  from 
superheated  steam  and  then  mixed  with  air  to  form  the  fuel/steam/air  reactants  mixture. 


Figure  36  ATR  reactor  with  steam  and  air  distribution  tubing 
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Figure  37  shows  the  detailed  design  of  the  ATR  reformer.  Because  catalyst  screening  was 
planned,  the  reformer  was  designed  for  ease  catalyst  replacement.  The  reformer  was  a  stainless 
steel  pipe  with  diameter  of  48.26  mm  and  length  of  47  cm.  The  reformer  was  placed  in  a  tubular 
furnace  in  order  to  keep  the  constant  operation  temperature  of  the  catalytic  zone.  Two  pieces 
stainless  steel  mesh  (mesh  size  20x20,  open  area  46%)  were  mounted  in  front  of  the  catalytic 
zone  aiming  to  distribute  the  mixing  flow  of  reactants.  Cordierite  monolithic  catalyst  was  placed 
in  the  catalytic  zone  and  its  length  was  60  mm.  After  reaction,  the  reformate  was  cooled  in  the 
copper  coils  as  shown  in  Fig.  35  and  the  residual  liquid  water  was  collected  in  a  condenser.  Then 
the  dry  reformate  was  analyzed  in  a  GC  system.  Agilent  GC6890  equipped  with  TCD  (thermal 
conductivity  detector)  was  employed  to  quantitatively  analyze  the  dry  reformate.  Molecular 
Sieve  5A  column  (Supelco,  60/80  mesh,  6FTxl/8  IN)  and  Hayesep  Q  column  (Agilent,  80/100 
mesh,  10FTxl/8  IN)  were  used  to  detect  molar  fractions  of  H2,  N2,  02,  CO  and  CO2.  The  total 
dry  reformate  volumetric  flow  rate  was  measured  by  Dwyer  flowmeter.  The  hydrogen  selectivity, 
carbon  selectivity  and  total  energy  conversion  efficiency  [136]  in  tests  were  defined  as  follows, 


Hydrogen  selectivity  (%)  = - - - xlOO 

(m  +  n/2)X)fy  fI 

(45) 

Carbon  selectivity  (%)=-^ — ^4-xlOO 

mxi|  H 

(46) 

LFTVT  +4,LHVro 

Efficiency  rj  (%)  =— - — - xlOO 

(47) 

For  n-dodecane,  in  the  above  equations  m  =  12  and  n  =  26;  and  for  desulfurized  Jet-A  fuel,  m  = 
1 1.3  and  n  =  22.6. 


Figure  37  Detailed  ATR  reactor  design 


4.5  Autothermal  reforming  experimental  tests  and  results 

In  the  ATR  tests  presented  in  section  4.5.1  to  4.5.4,  «-dodecane  (>99%)  was  used  as  the 
surrogate  fuel  of  Jet-A  to  study  the  autothermal  reforming  characteristics  of  the  reformer  and  the 
proposed  catalyst.  Rachner  [28]  reported  an  empirical  formula  for  Jet-A  to  be  C11.6H22.3  with 
molecular  weight  about  161.9  g/mol  and  LHV  of  43.26  MJ/kg.  As  a  surrogate  fuel,  /7-dodecane 
has  the  chemical  formula  of  C12H26  with  hydrogen  content  15.28  wt.%.  And  the  molecular 
weight  was  170.3  g/mol,  and  LHV  was  44.14  MJ/kg. 

4.5.1  Catalyst  screening  based  on  test  results  of  surrogate  fuel  (n-dodecane) 

The  major  advantage  of  using  monolith  as  the  catalyst  support  is  that  monolith  has  high  open 
frontal  area  so  that  the  flow  pressure  drop  is  low.  But  the  monolith  also  has  drawbacks  compared 
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to  fixed  bed.  The  mass  and  heat  transfer  rates  are  not  strong  due  to  low  Reynolds  number,  and 
the  mass  transfer  in  radial  direction  doesn’t  exist  in  monolith.  Therefore  the  high  loading  of 
catalyst  and  uniform  distribution  of  active  components  on  the  wall  of  the  monolith  are  crucial  for 
the  catalyst  preparation.  Six  catalyst  combinations  were  tested  and  compared  in  the  present 
study.  The  catalyst  samples  details  are  listed  in  Table  6.  Catalyst  1  to  3  had  higher  loading  of 
alumina  and  cerium  oxide  compared  to  the  rest.  Catalyst  4  to  6  have  different  loading  of 
potassium  oxide  and  nickel  oxide.  The  six  catalysts  were  tested  separately  in  the  ATR  reformer 
at  identical  operation  conditions:  flow  rate  of  «-dodecane  =  6  ml/min,  O2/C  =  0.45,  S/C  =  1.8, 
steam  was  pre-heated  to  210  °C,  air  was  pre-heated  to  165  °C,  «-dodecane  was  pre-heated  to  180 
°C,  and  the  operation  temperature  of  the  reformer  catalytic  zone  was  kept  at  700  °C.  Figure  38 
compares  the  mole  fraction  of  H2,  CO  and  CO2  in  dry  reformates.  Figure  39  indicates  the 
percentages  of  hydrogen  selectivity,  carbon  selectivity,  and  energy  conversion  efficiency. 

Figure  38  shows  that  catalyst  5  produced  the  highest  hydrogen  concentration  of  43.6%  as  well 
as  the  highest  carbon  monoxide  concentration  of  12.9%.  Since  both  hydrogen  and  carbon 
monoxide  are  fuel  sources  for  SOFC,  the  addition  of  hydrogen  and  carbon  monoxide  is  used  to 
compare  the  performance  of  different  catalysts.  Catalyst  5  had  the  highest  H2+CO  concentration 
of  56.5%,  followed  by  catalyst  4  with  52.8%  H2+CO  concentration  and  catalyst  6  with  48.6% 
H2+CO  concentration.  Figure  39  shows  the  same  tendency  of  hydrogen  selectivity,  carbon 
selectivity  and  energy  conversion  efficiency.  Catalyst  5  had  the  best  hydrogen  and  carbon 
selectivity  of  99%  and  88%  respectively.  The  efficiency  achieved  84.5%,  and  the  output  power 
was  2.8  kWt.  Catalyst  5  was  then  used  in  further  study  to  characterize  the  reforming  system. 


Table  6  Catalyst  screening  test  samples 

Catalyst 

number 

Composition  weight  ratio  (wt.  %) 

AI2O3 

Ce02 

La2C>3 

K20 

NiO 

Rh 

1 

7 

3 

2.5 

2.5 

5 

0.3 

2 

7 

2 

0.4 

0.8 

5 

0.3 

3 

7 

4 

0.4 

0.8 

5 

0.3 

4 

3.5 

1 

0.4 

0.8 

5 

0.3 

5 

3.5 

1 

0.4 

1.6 

5 

0.3 

6 

3.5 

1 

0.4 

0.8 

10 

0.3 

Figure  38  Mole  fraction  of  H2,  CO  and  CO2  in  dry  reformates 
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Catalyst  sample  number 

Figure  39  H2  selectivity,  carbon  selectivity  and  reaction  efficiency  comparison 

4.5.2  Effect  of  reformer  operating  temperature 

The  ATR  reaction  is  a  combination  of  exothermic  partial  oxidation  reaction  (POX)  occuring 
in  the  upstream  of  the  catalyst  sequentially  followed  by  endothermic  steam  reforming  reaction 
(SR).  Ideally  the  heat  released  by  POX  can  sustain  SR  if  the  reformer  wall  is  adiabatic,  and 
overall  the  ATR  reaction  is  thermally  neutral  or  slightly  exothermic.  However,  in  the  present 
study  the  ATR  reformer  was  not  thermally  insulated,  so  that  external  heat  source  was  needed  to 
sustain  the  catalytic  zone  at  constant  operation  temperature.  The  tubular  furnace  was  also  used  to 
heat  the  catalytic  zone  to  light-off  temperature  before  reactants  entering  the  reformer  since  POX 
reaction  does  not  occur  at  ambient  temperature. 

The  temperature  was  programmed  to  increase  from  400  °C  to  750  °C.  Seven  different 
temperature  points  were  tested  and  at  each  point  the  temperature  were  kept  for  40  min. 
Temperature  higher  than  750  °C  cannot  be  reached  due  to  the  endothermic  SR  reaction  and  the 
cooling  effect  of  the  reactants  flow  and  reformate  flow.  All  the  other  operation  conditions  were 
identical:  flow  rate  of  «-dodecane  =  6  ml/min,  O2/C  =  0.45,  S/C  =  1.8,  steam  was  pre-heated  to 
210  °C,  air  was  pre-heated  to  165  °C,  n-dodecane  was  pre-heated  to  180  °C.  Figure  40  shows  the 
product  molar  fractions  in  dry  reformate  at  different  temperatures,  and  Fig.  41  shows  the 
hydrogen  selectivity,  carbon  selectivity,  and  efficiency  at  different  temperatures.  It  can  be 
observed  from  Fig.  40  that  the  H2  molar  fraction  is  at  a  relative  stable  level  between  39%  and 
43.5%  at  different  temperatures,  whereas  CO  fraction  increases  from  5%  to  15%,  and  CO2 
decreases  from  15%  to  10%  as  the  temperature  increases  from  400  °C  to  750  °C.  CH4  has  a  high 
concentration  of  3.4%  at  400  oC  and  decreases  gradually  as  the  temperature  increases.  At  650  °C 
the  CH4  concentration  reduced  to  0.04%.  Hydrogen  selectivity,  carbon  selectivity,  and  efficiency 
are  all  quite  low  at  400  °C  as  shown  in  Fig.  41,  and  they  all  increase  gradually  as  the  temperature 
increases.  At  700  °C  the  hydrogen  selectivity  and  efficiency  are  slight  higher  than  that  at  750  °C, 
which  indicates  the  best  performance  of  the  ATR  reformer  can  be  obtained  at  700  °C  operation 
temperature. 


Project  final  report  August  5th,  2015 


400  500  600  700  800 

Operation  temperature  (  C) 


(a) 


Figure  40  (a)  H2,  CO,  CO2  and  CH4  fractions  in  dry  reformate  at  different  operation 

temperatures;  (b)  CH4  fraction  in  dry  reformate  at  Afferent  operation  temperature s 


Operation  temperature  (°C) 

Figure  41  H2  selectivity,  carbon  selectivity  and  reaction  efficiency  at  different  operation 
temperatures 
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4.5.3  Effects  of  S/C  (steam/carbon)  and  OfC  (oxygen  carbon)  ratios 

Steam  to  carbon  ratio  and  oxygen  to  carbon  ratio  are  crucial  variables  to  the  reformer  syngas 
compositions.  Qi  et  al.  [126]  reported  the  optimal  S/C  and  O2/C  for  ATR  were  in  the  range  of 
1. 5-2.0  and  0.35-0.5,  respectively.  By  keeping  other  parameters  as  constants  (flow  rate  of  n- 
dodecane  =  6  ml/min,  steam  was  pre-heated  to  210  °C,  air  was  pre-heated  to  165  °C,  n-dodecane 
was  pre-heated  to  1 80  °C,  and  the  operation  temperature  of  the  reformer  catalytic  zone  was  kept 
at  700  °C),  nine  different  S/C  and  O2/C  combinations  were  tested  in  the  reformer.  S/C  ratio  was 
controlled  by  water  flow  rate  and  O2/C  was  controlled  by  compressed  air  flow  rate.  Figure  42 
shows  the  H2  concentration  in  dry  reformates  at  various  S/C  and  O2/C.  Figure  43  shows  the  CO 
concentration  at  the  same  operation  conditions.  Both  H2  and  CO  concentrations  were 
significantly  influenced  by  S/C  and  O2/C.  Figure  42  shows  that  O2/C  has  relatively  stronger 
effect  on  H2  production  compared  with  S/C.  As  O2/C  increases  from  0.45  to  0.48,  H2 
concentration  has  a  tendency  to  decrease.  This  could  be  caused  by  the  combustion  of  FB  with 
excessive  O2.  As  seen  from  Fig.  43,  S/C  influences  CO  production  more  than  O2/C.  CO 
concentration  varies  at  a  range  of  1%  at  different  O2/C  conditions,  yet  it  can  change  at  a  range  of 
4%  at  different  S/C  operations.  The  potential  causes  could  be  combinations  of  CO  combustion 
and  water  gas  shift  reaction.  The  energy  conversion  efficiency  was  determined  by  both  the  H2 
and  CO  productions  and  the  total  production  flow  rate.  Figure  44  shows  the  efficiency 
comparison  at  different  operation  ratios.  It  can  be  concluded  that  the  efficiency  decreases  as  the 
S/C  increases,  and  at  S/C  =  1.5,  O2/C  =0.45  the  highest  efficiency  of  85.7%  is  achieved.  At  S/C 
=  1.5,  02  C  =0.42  and  O2/C  =0.48  the  efficiencies  of  85.5%  and  85.1%  are  slightly  lower  than 
the  highest  value. 


Steam  to  carbon  ratio 

Figure  42  H2  concentration  in  dry  reformate  at  various  S/C  and  O2/C 
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Staam  to  carbon  ratio 

Figure  43  CO  concentration  in  dry  reformate  at  various  S/C  and  O2/C 


Steam  to  carbon  ratio 

Figure  44  Energy  conversion  efficiency  at  various  S/C  and  O/C 

4.5.4  Coke  formation 

Figure  45  compares  the  carbon  selectivity  of  CO  and  CO2  at  various  S/C  and  O2/C.  It  is 
noticed  that  in  the  high  energy  conversion  efficiency  operation  conditions,  the  carbon  selectivity 
is  relatively  low  in  the  range  of  82%  to  87%.  Since  the  detectable  CH4  wTas  below  0.5%  for  all 
tests,  it  was  believed  there  was  high  possibility  that  the  unconverted  carbon  in  the  fuel  was 
transferred  to  coke.  After  about  80  hours  running  of  the  ATR  system,  the  fuel  heating  hose  was 
dissembled  and  it  was  found  that  coke  formation  in  the  tube  fittings  as  well  as  inside  the  heating 
hose  was  quite  severe  at  both  ends  of  the  heating  hose.  Figure  46  shows  pictures  of  coke 
formation  in  the  two  fittings.  The  fitting  at  the  front  end  had  more  serious  coke  problem  than  the 
one  at  the  rear  end.  The  downstream  path  of  the  heating  hose  was  tested  with  N2  flow  and  no 
apparent  coking  was  detected.  The  front  cap  of  the  ATR  reformer  was  dissembled  to  check  for 
coke  formation  and  both  the  flange  cap  and  fuel  pressurized  nozzle  were  free  of  coke  as  shown 
in  Fig.  47.  It  draws  the  conclusion  that  the  most  severe  coke  formation  was  caused  by  the  high 
pre-heat  temperature  for  the  liquid  fuel. 
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Figure  45  Carbon  selectivity  at  various  S/C  and  OFC 


fitting  at  front  end  fitting  at  rear  end 


Figure  46  Coke  formation  in  heating  hose  fillings 


Figure  47  Flange  cap  and  nozzle  are  free  of  coke 

In  order  to  control  the  coke  formation  in  the  ATR  system,  the  fuel  pre-heat  temperature  was 
reduced  and  the  pre-heat  temperatures  of  steam  and  air  were  increased  to  supply  more  heat  for 
fuel  evaporation.  Six  tests  of  different  pre-heating  temperatures  were  conducted  and  the  details 
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are  listed  in  Table  7.  Test  #7  is  the  test  with  unacceptable  level  of  coke.  Temperatures  higher 
than  250  °C  for  steam  and  175  °C  for  air  were  difficult  to  achieve  because  of  limitation  of  the 
heating  elements  and  significant  heat  loss  along  the  transportation  paths.  Other  operation 
conditions  were  identical  and  at  the  optimized  conditions:  flow  rate  of  n-dodecane  =  6  ml/min, 
O2/C  =  0.45,  S/C  =  1.5,  and  the  reformer  temperature  =  700  °C.  Details  of  the  temperatures  tests 
as  well  as  the  mole  fraction  of  substances  in  the  dry  reformates  are  listed  in  Table  7.  Figure  48 
shows  the  energy  conversion  efficiency  and  carbon  selectivity  at  different  temperature 
conditions.  It  can  be  seen  that  the  H2  molar  fraction  and  the  energy  conversion  efficiency 
increased  as  the  pre-heating  temperatures  of  reactants  increased.  When  the  fuel  temperature 
reached  100  °C,  the  highest  carbon  selectivity  of  98%  was  achieved  and  the  corresponding 
efficiency  was  79.7%.  As  the  fuel  temperature  increased  to  120  and  140  °C,  the  carbon 
selectivity  reduced  to  96%  and  93%  respectively,  but  the  energy  conversion  efficiency  increased 
to  81.7%  and  83.5%.  Considering  both  the  carbon  selectivity’  and  energy  conversion  efficiency, 
pre-heating  temperatures  in  test  #6  was  selected  as  the  best  operation  conditions.  After  20  hours 
running  at  this  condition,  the  fuel  heating  hose  was  dissembled  from  the  ATR  system  and  no 
significant  coke  formation  was  observed  at  the  fittings  in  both  ends. 


Table  7  Tests  of different  pre-heat  temperature^ 


Fuel 

Steam 

Air 

Mole  fraction  (%) 

Test# 

temperature 

(°C) 

temperature 

(°C) 

temperature 

(°C) 

h2 

CO 

C02 

ch4 

1 

80 

250 

155 

33.70 

15.80 

9.70 

0.90 

2 

80 

250 

165 

35.40 

15.98 

8.90 

0.61 

3 

80 

250 

175 

35.79 

15.48 

10.45 

0.64 

4 

100 

250 

175 

35.80 

15.60 

10.83 

0.41 

5 

120 

250 

175 

37.48 

15.30 

10.13 

0.64 

6 

140 

250 

175 

38.66 

15.35 

9.15 

0.42 

7 

180 

210 

165 

42.80 

13.94 

10.11 

0.12 

Figure  48  Efficiency  and  carbon  selectivity  of  temperature  tests 
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4.5.5  Reforming  of  desulfurized  Jet- A  fuel 

Sulfur  clean  real  commercial  Jet-A  fuel  was  tested  at  the  optimized  and  coke  formation 
reduced  ATR  conditions.  The  sulfur  concentration  in  the  tested  desulfurization  Jet-A  fuel  was 
about  20  ppmw.  Figure  49  shows  the  mole  fraction  of  H2,  CO,  CO2  and  CH4  in  the  dry  reformate 
at  reforming  time  of  0,  30,  60  and  90  min.  It  can  be  seen  that  the  production  of  H2  and  CO  is 
relatively  stable.  The  mole  fraction  of  H2  is  around  35%  and  the  mole  fraction  of  CO  is  about 
12.5%.  Concentration  of  CO2  is  in  the  similar  level  with  CO.  ATR  of  desulfurized  commercial 
Jet-A  fuel  produces  more  CH4  than  ATR  of  n-dodecane.  Averagely  2%  CH4  is  produced  in  ATR 
of  real  Jet-A  fuel  while  it  is  usually  below  1%  in  ATR  of  «-dodecane.  Figure  50  shows  the 
average  energy  conversion  efficiency,  H2  selectivity  and  COx  selectivity  of  ATR  of  desulfurized 
commercial  Jet-A  fuel.  The  H2  selectivity  is  93%  and  the  COx  selectivity  is  above  96%,  so  that 
the  coke  formation  should  not  be  a  severe  problem.  However,  the  energy  conversion  efficiency 
for  the  ATR  of  desulfurized  commercial  Jet-A  fuel  is  lower  than  that  for  the  ATR  of  rc-dodecane 
at  the  same  operation  conditions.  For  rc-dodecane  the  energy  conversion  efficiency  is  83.5% 
while  for  desulfurized  commercial  Jet-A  fuel  the  efficiency  reduces  to  77%  with  all  identical 
operation  conditions. 


Figure  49 


(a)  Mole  fraction  of  products  in  the  dry  reformate;  (b)  Mole  fraction  of  CH4  in  the 

dry  reformate 
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Figure  50  Energy  conversion  efficiency,  H2  selectivity  and  COx  selectivity 

5.  Fuel  adaptability  study  of  a  lab-scale  2.5  kWth  autothermal  reformer 

Fuel  adaptability  of  a  lab-scale  2.5  kWth  autothermal  reforming  system  was  experimentally 
investigated.  Analyses  of  thermodynamics  of  autothermal  reforming  of  different  fuels  at  various 
temperatures  were  provided  to  determine  the  associated  flow  rates  of  fuel,  air,  and  steam.  The 
tested  reformer  has  a  specifically  designed  mixing  chamber  and  NiO-Rh  based  bimetallic 
catalyst  with  promoters  of  Ce,  K  and  La,  which  were  originally  designed  for  autothermal 
reforming  of  Jet-A  fuel.  The  adaptability  of  the  reformer  to  multiple  fuels,  including  diesel, 
gasoline,  ethanol,  and  methanol  was  experimentally  studied  at  different  operating  conditions. 
Molar  fractions  of  H2,  CO,  CO2,  and  CH4  in  the  dried  reformate  were  measured  and  the  results 
were  compared.  The  autothermal  reforming  performance  of  the  reformer  for  different  fuels  was 
evaluated  with  regard  to  hydrogen  selectivity,  COx  selectivity,  and  energy  conversion  efficiency. 
It  was  concluded  that  the  currently  developed  reformer  and  the  associated  catalyst  have  great  fuel 
adaptability  and  high  potential  being  used  in  a  fuel-flexible  autothermal  reforming  system. 

Logistic  fuels  of  gasoline,  diesel,  and  Jet-A  fuel,  and  potential  alternative  fuels  of  methanol 
and  ethanol  were  tested  in  the  lab-scale  ATR  system.  Gasoline  with  total  sulfur  concentration  of 
9.6  ppmw,  and  diesel  with  total  sulfur  concentration  of  5.4  ppmw  were  purchased  from  local 
Shell-Oil  gas  station.  Jet-A  fuel  with  initial  total  sulfur  concentration  of  1140  ppmw  was 
purchased  from  local  airport  fuel  service  and  was  desulfurized  to  20  ppmw  by  homemade  NiO- 
Ce02/Al203-Si02  adsorbent  prior  to  reforming.  Important  properties  of  the  fuels  in  the  reforming 
tests  are  listed  in  Table  8. 

In  Fig.  51  to  Fig.  55  the  molar  fractions  of  reformate  for  five  different  types  of  fuels  are 
presented.  The  reforming  operating  conditions  including  inflow  molar  ratios  of  O2/C,  S/C 
(steam/carbon),  total  flow  rate,  reforming  temperature,  and  GHSV  (gas  hourly  space  velocity) 
are  given  for  each  tested  case. 


Table  8  Properties  of  the  tested  fuels  (Jet-A  desulfurized  in-house) 


Gasoline  Diesel  Jet-A  Methanol  Ethanol 


Average  chemical  formula 

CaHix 

Ct*Hj6 

CiLtHara 

CHiOH 

CjHsOH 

Molecular  weight  (g  mol-1) 

114 

194 

161.5 

32 

46 

Boiling  point  © 1  atm  (  C) 

25  205 

180  360 

176 

66 

79 

liquid  density  C  IS  C  (kg  m~3) 

719.7 

832 

800 

791.8 

789 

Lower  heating  value  (MJ/kg) 

44.4 

43.4 

43.26 

19.93 

28.86 

Sulfur  content  (ppmw) 

9.6 

5.4 

20 

0 

0 
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Figure  51  Results  of  AT R  of  gasoline 
(OfC^O.4,  S/C=2,  reformer  temperature 
= 700  °C,  fuel  flow  rate=6.1  ml/min, 
GHSV=  41,370  h'1). 


Figure  53  Results  of  AT R  of  Jet-  A 
(02/0=0.45,  S/C=1.5,  reformer  temperature 
=700  °C,  juel  flow  rate=6.0  ml/min, 

GHSV=  35,980  Kl). 


35 


5 


o  y  .  i  — i - . - 

0  20  40  60  60  100 

Reforming  time  (min) 

Figure  52  Results  of  AT R  of  diesel 
(O2/C=0.4,  S/C=2,  reformer  temperature 
= 700  °C,  fuel  flow  rate  =5. 4  ml/min, 
GHSV=  42,850  K1). 


Figure  54  Results  of  AT R  of  methanol 
(O2/C=0.2,  S/C=1.7,  reformer  temperature 
=  700  °C,  fuel  flow  rate=12.4  ml/min, 
GHSV=59.250  K1). 


Figure  55  Results  of  AT R  of  ethanol  Figure  56  Fractions  of  FI 2,  CO,  CO2,  and 

(02/0=0.35,  S/C=3,  reformer  temperature  CH4  in  the  dry  reformate. 

=700  °C,  fuel  flow  rate=8.6  ml/min, 

GHSV=62,630K‘). 
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Figure  57  Comparison  of  FI 2  selectivity,  COx 
selectivity,  and  energy  conversion  efficiency 


Figures  56  and  57  show  the  comparison  of  K2  selectivity,  COx  selectivity,  and  energy 
conversion  efficiency  of  five  of  the  tested  fuels.  Since  the  reformer  is  particularly  developed  for 
Jet-A  fuel,  better  reforming  performance  is  seen  for  Jet-A.  However,  in  general,  the  reformer  can 
also  be  adaptable  to  other  four  types  of  fuels. 


6.  Integrated  system  of  autothermal  reforming  and  shifting  with  heat  recovery 

The  integration  of  a  hydrocarbon  fuel  auto-thermal  reformer  with  a  water- shift  reactor  as  one 
compact  reactor  has  been  designed  and  fabricated,  as  schematically  shown  in  Fig.  58.  The 
reactor  include  three  sections — reforming,  heat  exchange  for  giving  heat  for  fuel,  air,  and  steam, 
and  water  shifting.  The  entire  unit  allows  accomplishment  of  thermal  management  to  have 
favorable  high  temperatures  for  autothermal  reforming  and  low  temperatures  for  water-shift 
reaction,  while  heat  from  reformate  is  recovered.  In  overall,  high  hydrogen  conversion  rate  from 
hydrocarbon  fuels  is  achievable  through  the  reactor. 


Reforming 


fu?4 


water 


ter  [ 
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- J — ~ — 1 

water  air  fuel  water 

4  111 


Shifting 


h2 


Figure  58  Schematic  of  integrated  unit  of  reformer  and  water-gas-shift 

The  auto-thermal  reforming  works  at  temperatures  from  500  °C  to  700  °C.  The  reformed  gas 
temperature  is  then  lowered  to  350  °C  in  four  sections  of  heat  exchange,  which  at  the  same  time 
heats  up  the  fuel,  air,  and  water  going  into  the  front  of  the  reformer.  The  water  shift  section 
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performs  at  temperatures  of  around  350  °C.  Test  results  are  expected  to  be  available  in  two 
months.  Supplemental  report  will  then  be  provide  to  UTC  and  NOR. 
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